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SUMMARY 


This report Ls concerned with development and extension of a computer 
program for the study of mission performance of spacecraft using solar electric 
propulsion. The code is designed for analysis of the NASA /Lewis SERT-C design 
which involves attitude constraints, however, it has been generalized to consider 
any of the following configurations: (1) yaw motion only (SERT-C), (2) yaw and 
roll only, (3) unconstrained motion. 

The code calculates time optimal or nearly time optimal geocentric tra- 
jectories. It uses the method of averaging to reduce computer time by orders of 
magnitude compared to a precision integrated trajectory. The averaged rates of 
change of the mean values of the state and costate are found by numerical quadrature. 
The differential equations for the mean state and costate may then be integrated in 
large time steps (typically days). A set of nonsingular orbital elements is used to 
avoid numerical difficulties for eccentricities and inclinations at or near zero. 

The radiation model is analytic to reduce run time. Included in the code are 
options for consideration of oblateness, solar motion, shadowing with or without 
a delay in thruster startup, an analytic radiation and power degradation model, and 
an initial high thrust stage of one or two impulses of specified ^V. 

A costate formulation of the problem yields a two point boundary value 
problem which is solved by a Newton iteration on the initial costate and value of 
transfer time. Initial values of the unspecified states and costates and a guess 
for the transfer time are chosen. An optimal trajectory is then generated by 
integrating the averaged state and cpstate using a Runge-Kutta method. A Gaussian 
quadrature averages the state and costate derivative over an orbit. This will 
generate an optimal trajectory to the wrong terminal state. A sensitivity matrix is 
then generated by varying the initial conditions and running a set of neighboring 
trajectories. A Newton iteration on the initial conditions is then used to drive the 
terminal errors to within specified bounds. The final converged trajectory is a 
minimum time trajectory (nearly minimum time for attitude constraints) for the 
specified velocity increment in the high thrust phase (if included). 
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The attitude constraint causes power to become a function of thrust direction 
and sun direction, and the time optimal thrust direction becomes a complex function 
of primer vector direction. The analysis for the attitude constrained case is 
considerably more complicated than for the unconstrained case. 

For the zero roll and pitch case a suboptimal control is developed which is 
nearly time optimal, uses less fuel than the time optimal, and is obtained from the 
solution of a cubic equation. The solution may yield discontinuous changes in the 
thrust direction. The time delay in thruster startup is modeled as a quadratic 
function of time in shadow. The proton and electron radiation field is modeled as 
equivalent 1 MEV particle flux as a function of spacecraft altitude and latitude as 
well as the solar cell shield thickness. The power loss function is modeled as an 
analytic function of fluence, cell thickness and base resistivity. Spacecraft 
parameters such as yaw angle and sun incidence angles may be displayed. 

Numerical results were obtained for many SERT-C and other type missions. 
Attitude constraints increase flight times by a few percent for SERT-C type missions. 
Power degradation can be quite severe at lower altitudes nearly doubling transfer 
time, compared to a no degradation case. 
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SECTION 1 


INTRODUCTION 


The use of time optimal trajectories can save considerable flight time in 
low thrust missions (days, weeks or even months). As a result there has been 
considerable use of optimization theory in the analysis of low thrust trajectories 
and missions. Recent advance in the state of the art have made possible routine 
calculation of geocentric optimal trajectories considering more representative 
environmental constraints and influences. It has recently been proposed^^^ that 
optimal trajectory computer programs form the basis of a guidance technique. 

Most of the geocentric trajectory analyses for solar electric spacecraft 

have assumed fully articulated spacecraft where the solar array and the thrust 

(2-5) 

vector can both be pointed in their respective optimum direction . The 
(fi ) 

SERT-C spacecraft study performed at NASA-Lewis introduced attitude motion 
constraints whereby the spacecraft was permitted to rotate only about an axis 
parallel to the Earth radius vector. This design constrains the thrust vector to 
lie in a plane perpendicular to the radius vector. The solar panels can rotate about 
an axis that is orthogonal to both the radius vector and the nominal thrust vector. 
Although the constraints simplified the SERT-C design and reduced the attitude 
sensor requirements^ , the constraint couples the power developed to the thrust 
direction. Most thrust orientations do not allow maximum power to be developed. 
The flight time is increased both because of this coupling and because of the thrust 
direction constraint. 

The present study, performed for the NASA Lewis Research Center, was 
undertaken both to provide software for the evaluation of the SERT-C design with 
attitude constraints and to provide a basis for the proposed guidance scheme. The 
effort described in this report builds on the work of References 2, 7, and 9. For 
that effort a computer program was developed which calculated minimum time 
trajectories for unconstrained solar electric and nuclear electric propulsion. Also 
an initial high thrust stage of specified AV with one or two impulses and a final 
impulse could be included. The low thrust stage uses Kryloff- Bogoliuboff 
averaging^^^^ of both the state and the costate. The averaged rates of change of 
the mean values of the state and costate are found by numerical quadrature. The 
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differential equations for the mean state and costate may then be integrated in 

large time steps (t 3 rpicaily days). A set of nonsingular orbital elements, the 

equinoctial elements^^ ^ is used to avoid numerical difficulties. The method of 

(12 13) 

averaging has been used extensively in recent years. Edelbaum " has used 
averaging to calculate analytic solutions for special cases of optimal low thrust 
trajectories, and others have used averaging when considering effects such as 
oblateness, third body perturbations and non-optimal thrusting^^"^’ Jasper^^^ 

utilizes equinoctial orbital elements and averaging in recent low thrust optimiza- 
tion work. An interesting feature of his work is numerical comparisons of 
averaging with integration of the full state and costate. He does not include the 
effects of perturbing forces, other than thrust, on the costate. 


The first stage high thrust optimization is based on a very efficient computer 
program developed by Huntington SmalP ' . This program uses a special set 

of variables and form of the switching conditions developed by Small. The initial 
orbit is assumed circular with specified semimajor axis and inclination, while the 
final orbit has specified semimajor axis, eccentricity, and inclination. During the 
high thrust phase an inverse square gravitational field is assumed. Because the 
initial orbit is circular, it was possible to constrain the initial costate to the region 
that yields solutions. This program rapidly calculates either one or two impulse, 
minimum-fuel, time-open trajectories. Because this transfer always requires 
less than a full revolution, its time is negligible compared to the low thrust phase 
and is not considered. A summary of the high thrust analysis is given in Appendix A. 


The utilization of high and low thrust in combination has been considered by a 

, . (19, 20, 21, 22, 23) 

number of authors 


The effect of oblateness is included by analytically adding its associated 
rate of change of the mean state and costate to that due to thrust. The effects of 
shadowing are calculated by assuming that thrust is turned off in shadow. The 
shadow entrance and exit times are calculated analytically by solving a quartic 
equation. The effects of radiation degradation are calculated by modeling equiva- 
lent 1 MEV electron flux as a function of radius and latitude. The power is then 
expressed as a function of the total accumulated particle fluence. The model in the 
first version of the program was valid for one cell and shield thickness, A new 
model, developed herein, is valid for a variety of cell thickness, shield thicknesses 
and base resistivities. As for ail perturbations, the effect of radiation degradation 
on the costate as well as the state is calculated. 


The overall trajectory is optimized by a shooting method. Initial values of 
the unspecified states and costates are chosen at the initial time. An optimum high 
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and low thrust trajectory is then generated by integrating the state and costate 
through both stages. This will generate an optimal trajectory to the wrong 
terminal state. A sensitivity matrix is then generated by varying the initial condi- 
tions and running a set of neighboring trajectories. A Newton iteration on the 
initial conditions is then used to drive the terminal errors to within specified 
bounds. The final converged trajectory is a minimum time trajectory (nearly 
minimum time for attitude constraints) for the specified velocity increment in the 
high thrust phase. 

The principle modifications and additions performed under this contract to 
the earlier work are: nearly time optimal trajectories can be generated for the 
attitude constrained case of zero roll and pitch or zero pitch and free roll as well 
as the unconstrained case; the shadow model has been modified to include the effect 
of a delay in thruster turn-on after leaving the shadow; a new, more accurate and 
more general solar array degradation model is developed; and additional spacecraft 
parameters may be calculated and displayed. All of the options of the earlier code 
are included except that the final impulse option has been removed. 

The problem, then, which is considered in this report is the calculation of 
time optimal or nearly optimal geocentric transfers using solar electric space- 
craft which may have attitude constraints and an optional initial high thrust stage. 
For the low thrust phase the initial and final orbits are general ellipses. An inverse 
square gravity field with oblateness is assumed. Thrust is assumed proportional 
to power with constant specific impulse and efficiency. The effect on power of 
solar distance may be included. The thrusters may be assumed to be off while the 
spacecraft is in Earttfs shadow and for a start-up delay time which corresponds to 
the sum of the time for the solar array to achieve operating temperature and the 
time for the thruster to achieve full thrust after the solar array power is applied 
to the power processor. The Van Allen radiation is modeled analytically and its 
effect on power degradation included. 

For the attitude constrained case the class of spacecraft modeled is indicated 
in Figure 1-1 where the roll axis lies in the orbit plane and is perpendicular to the 
Earth-spacecraft radius vector; the yaw axis is parallel to the Earth radius vector 
directed toward the Earth; the pitch axis is perpendicular to the orbit plane and 
directed south. For the nominal attitude the principle body-centered axes are 
aligned with this coordinate system. The ion thrusters are mounted on the negative 
roll face of the spacecraft and directed parallel to the roll axis. The solar arrays 
are flat panels and are capable of rotation about their longitudinal axis, which is 
aligned with the spacecraft pitch axis. The required low thrust directions are 
achieved entirely by the spacecraft attitude rotations. 
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Figure 1- 1 Spacecraft Configuration 


Ttie thrust vector must have no component along the vector from Earth to 
spacecraft. The thrust direction is determined by one control variable, the yaw 
angle. Since pitch and roll are constrained to be zero, the orientation of the axis 
of the solar panels is determined by the yaw angle. The panels are allowed to 
rotate about this axis, resulting in one more control variable. Power is assumed to 
be proportional to the cosine of the angle between the normal to the panels and the 
line to the sun. The initial state is given and it is desired to reach a final subset 
of the state at some, unspecified, final time. Although the original aim was to 
derive the time optimal control, instead, a particular suboptimal control, which is 
nearly time optimal, is derived. A costate formulation is used along with the method 
of averaging to set up a two point boundary value problem which can be solved by 
using a Newton method to iterate on the unknown initial costate and value of the final 
time in order to meet the desired final conditions on the state and costate. 

The computer program developed is used to assess the increase in flight 

time and fuel consumption due to introducing the attitude constraints and other 

perturbations. A User's Manual which contains listings and descriptions of the 

(24) 

subroutines and a description of the use of the code is published separately 

(25 26) 

Two papers have been presented based on the material in this report ' 

Reference 25 has been accepted in modified form to be published in The Journal of 
Spacecraft and Rockets. 
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SECTION 2 


GENERAL TECHNIQUE 


2. 1 Introduction 

Three areas of general technique are discussed In this section. One Is the 
method of averaging, essential to the running of a program which generates many 
trajectories in a reasonable amount of computer time. Second is discussed the 
method of generation of the final trajectory. A time optimal or nearly optimal 
trajectory is desired, A state and costate formulation is used which results In a 
two point boundary value problem which can be solved by a Newton iteration 
procedure. Finally some comments on numerical techniques are made. Later 
sections give details of the dynamical equations, calculation of the control, and the 
effects of perturbations on an individual trajectory. 


2. 2 Averaging 

A great savings in computer time can be effected by considering a first 
approximation to the state and costate. Short period variations In the state and 
costate are eliminated by the averaging technique. When low thrust propulsion Is 
utilized and the other perturbations to the inverse square motion are small and 
when the state Includes the five slowly varying orbital elements which indicate the 
size, shape and orientation of an orbit and possibly other slowly varying quantities, 
then averaging may be used. Spacecraft mass and the accumulated particle fluence 
(both slowly varying) were also averaged. The orbital element indicating the 
position of the spacecraft In the orbit Is eliminated by the averaging process. 

The averaged Hamiltonian can be defined as 


iS 



H dt 


( 2 . 1 ) 


where H is the unaveraged Hamiltonian and T is the orbital period. When calcu- 
lating this Integral the state and costate are held fixed. The motion of the space- 
craft is assumed to vary in a manner described by Kepler's equation over the 
averaging integration. The approximate state and costate satisfy the canonical 
equations 
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( 2 . 2 ) 



(2.3) 


where the overbar LndLcates the approximate quantities. 

In what follows the averaging integral for oblateness (J^) is solved analyti- 
cally; otherwise a numerical quadrature formula is used. The differential equations 
can then be solved numerically using a time step which is much larger than but 
unrelated to the number of orbital revolutions. 


2. 3 Optimal and Suboptimal Trajectories 

The original aim was to calculate time optimal trajectories for the constrained 
problem as well as the unconstrained problem. However in the course of the 
analysis, which is described more fully in Section 3, it was decided to use a slightly 
suboptimal control for the constrained case of zero roll and pitch. The suboptimal 
control chosen is nearly time optimal, saves on fuel, and it is easier to calculate. 

The controls considered for the zero pitch but free roll case and the unconstrained 
case are time optimal. Transversality conditions for the time optimal control are 
used for all cases. In the following discussion the words ’^optimal trajectory" will 
be used to refer to both the optimal and suboptimal trajectories. In all cases the 
method used to generate a low thrust trajectory is to develop the Hamiltonian, to 
calculate a control (the thrust direction), either optimal or suboptimal, and to write 
the canonical equations for the state and costate. Some or all of the initial state 
and costate are specified. Application of transversality conditions for the time 
optimal problem yield additional specifications on the state and costate. Thus a two 
point boundary value problem results which must be solved to obtain the requisite 
trajectory. When these equations are solved, an extremal trajectory will result 
which is usually locally optimal. No attempt is made to investigate generalized 
Jacobi-type conditions to establish local sufficiency. Also, in common with other 
nonlinear problems, there may be more than one extremal meeting the same end 
conditions. The very difficult question of global optimality is not considered. 

The single trajectory generation portion of the code is coupled with a 
Newton iterator to solve the two point boundary value problem. The unknown initial 
conditions and value of the final time are iterated on in order to meet the final 
conditions which are functions of the final state and costate (usually a specified 
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semLmajor axis, eccentricity, and inclination and the appropriate transversality 
conditions). The partial derivative matrix of final conditions with respect to the 
initial costate is obtained numerically by calculating neighboring trajectories to 
a nominal. 

When only low thrust is included the initial orbit is specified as well as the 
mass of the spacecraft. A seventh state variable, the equivalent 1 MEV electron 
flug];^ce is also specified. At the unknown final time either five orbital elements 
can be specified or three: the semimajor axis, a, the eccentricity, e, and the 
inclination, i. Transversality conditions require that the adjoint to the mass and 
fluence be zero at the final time. If the final line of nodes (fi) and argument of 
perigee (aj) are free, then their adjoints are zero (X^ " either 

case the final value of the Hamiltonian must be unity. 

If initial high thrust is included, the necessary conditions for minimum time 
low-thrust transfer with specified high thrust velocity increments can be derived 
by considering each phase separately as a variational problem. The cost of each 
phase can then be expressed as a function of its terminal states and costates. The 
proper interface conditions for each phase can then be derived by considering tne 
parameter optimization problem of minimizing the time of the low thrust phase for 
fixed velocity increments in the high thrust phases. This minimization is carried 
out over all the free states and costates at the interfaces between high and low thrust. 

The flight time for the high thrust stages is assumed to be negligible 
compared to the low thrust stage flight time. At the beginning of the first high 
thrust stage a, e, i are assumed given (eccentricity is set to zero). Transversality 
conditions then require Xq = 0 and 

For the high thrust phases only an inverse square gravity field is assumed. 

X is a constant of the motion and therefore X^ remains zero during the initial 
hl^ thrust phase if Q is fixed. A single extremal is generated in the following 
fashion (it does not necessarily reach the desired final conditions). Values for Xg^. 

X , X-. and ware picked. Values for X^^. X^. are also picked but not used until 
after the first high thrust stage. Since the initial orbit is circular we are assuming 
that the first impulse occurs at u = f + wwhere f = 0 and sowjust indicates the loca- 
tion in the orbit of the first impulse. From the standpoint of maximizing the primer 
vector magnitude, the location on tne orbit at wh ich an impulse occurs is arbitrary, 
since for a circular orbit one position cannot be distinguished from another. 
Optionally, the user may specify the initial line of nodes. There is then an impulse 
in the direction of the primer vector. The magnitude of the ^V will be equal to the 
maximum for an optimal one impulse transfer if that maximum is less than the 
specified total initial ^V, or otherwise equal to the specified total initial ^V . If 
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the maximum one Impulse is less than the total available, then a second 
impulse occurs at the next maximum of the Hamiltonian, again in the direction of the 
primer vector. The program requires that all remaining initial high thrust be 
used for this impulse. This may not be optimal, as one or more additional impulses 
may be optimal. Generally, for practical cases one or two impulses will be optimal. 

After the one or two initial impulses there is a resultant orbit a, e, i, f) , 
a' and values for X„» X:» X,-.» X ■ The high thrust program scales the costate 

& vJ L U) 

so that the magnitude of the primer vector is 1 at the impulses. At the beginning 
of the low thrust stage they may need to be rescaled so that final transversality 
conditions will be met. This scaling factor, the picked values of x^, Xj^ and the 
values of X„, X, X., X_ , X from the high thrust stage are then used as input to 
the low thrust equations of motion. These are essentially Eqs. (2. 2) and (2. 3). 

They are integrated to the picked final time, t^. The final orbit is then reached. 

This is an extremal trajectory, though possibly not arriving at the desired final 
orbit. 


We desire to reach a specified a, e, i. Transversality conditions require 
that , X^, X^, Xj^ be zero and that H = 1. If we have not met these final con- 
ditions a Newton iteration on the free initial conditions and time of flight is used in 
order to meet the final conditions to within some tolerance. 


The actual calculations of the initial high thrust stage uses code developed 
(12 13) 

by H. Small ' which uses a special set of variables. (See Appendix A. ) The 
low thrust code uses equinoctial orbital elements rather than classical orbital 
elements. This stage is discussed in detail in Section 3. 


The Newton method works by first guessing values for the iteration parameters, 
call them x and t^, and then running a nominal trajectory which will yield final 
conditions y which in general are not equal to the desired final conditions, y^. 

Revised values for x, t^ may then be obtained by calculating a sensitivity matrix or 
partial derivative matrix. A, which is generated by varying slightly, one at a time, 
each of the iteration parameters, x, and running a new, neighboring, trajectory. 
Differencing the resulting values of the final conditions with the nominal values 
yields a ^y for each In addition can be calculated analytically except perhaps 
for which can be approximated numerically by varying t^ slightly and evaluating 
the corresponding H, differencing this with the nominal H and dividing by ^t^. Then 
A is an approximation for the partial of y with respect to x, t^. 


A 




(2.4) 
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A revised estimate of the iteration parameters can then be obtained by the formula 


= f] - A'\y-y^) (2.5) 

NEW L ^JoLD 

A new nominal trajectory can then be generated and the procedure continued until 
the final conditions are met to within some tolerance. In the event that the new 
X, t^ do not yield a reduction in the norm of the final condition errors, the change 
in (x, t^) is reduced in magnitude by factors of 2. Also there is an option of using 
a modified Newton-Raphson procedure wherein the A matrix is not always recalcu- 
lated at each iteration by running neighboring trajectories, but instead a new A 
may be approximated using the old A and the values of the changes in x, t^ . 

We will now specify more precisely the initial and final conditions. Since the 
low thrust code uses equinoctial elements, they will be defined here in terms of 
the classical elements, 

a = a 

h = esin(4o+Q) 

k = ecosC^ij+Q) (2.6) 

p = tan^-^^sinp) 

q = tan ^ cos Q. 

These are five of the seven state variables. The other two are mass, m, and 
particle radiation fluence, N. 

The Newton iteration scheme iterates on a vector x to drive a vector y to 
zero, where now we are using y to mean y - y^ as used above. Separately coded is 
the iteration variable t^, the final time, x is a function of the unknown state or co- 
state elements at the initial time; y is the error in tlie final conditions, a function of 
the state and costate at the final time and the Hamiltonian at the final time. 

Two final conditions options are considered. In the first, all five orbital 
elements are specified at the final time. Thus, 

a(tj) - 
h(tf) - 

k(t^) - 
P(tf) - 
q(tf) - 

- 1 





P 



where the subscript ^ indicates the desired value. 

For the second option only three orbital elements (a, e, i) are specified. 
The adjoints to n (t^) and must then be driven to zero. They are given in 

terms of the equinoctial state and costate. 

= X + qx„ - px„ (2. 8) 

n w p q 

X = kx. - hx^ ( 2 . 9) 

CO « k 

The eccentricity and inclinations are given by 


e = Jh^ - 

L = 2tan'^yp^ + q^ 

For this case, then, let 


'^p^(tj) + q^(tj.) - tan-^ 

h(tj)Xj^(tf) - k(tj)Xj^(tj) 
P(tf)Xq(tf) - q(t^)Xp(tj) 

H(t^) - 1 


a(tf) - a^ 

*^h^(t^) - k^(tj) - e^ 


( 2 . 10 ) 
( 2 . 11 ) 


{ 2 . 12) 


When 0 driving the above y to zero is equivalent to driving p and q to zero and 
not considering pX^ “ k\p. When calculating small errors near above y 

expression cannot distinguish between positive or negative inclinations. This can 
lead to numerical problems, so when i^ = 0, the code sets the third component of 
y to p(t^) and the fifth component to q(t^) in the SEP code. A similar situation exists 
if e^ = 0, so in this case the second component is h(t^) and the fourth component 
k(t^). If both ” 0 and « 0, then both options are equivalent in the SEP code. 
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For low thrust only, the Initial state must be specified, and therefore x is 
just the initial costate. When initial impulses are included either and q are 
free or just ^ is free. Since the initial orbit is circular for this case, indicates 
the location of the first impulse. Thus \ (t ) = 0 . If o (t ) is free, \ (t ) = 0. 

60 O O ^ 

It is shown in Appendix A that for the Q (t^) specified case, setting the initial 

determines the initial \ , Thus is used as the iteration variable. The high 

thrust code requires special input variables and these variables must be within 
specific bounds if a circular parking orbit is assumed (Appendix A). Taking into 
account these bounds let f ^ first three iteration parameters where 

Small’s special variables are given by 


T 



(2. 13) 


k = COST (. 75 -f , 25 




^2 


j = (1+k COS T ) J 


COS T‘k 

cos T+k 




(2. 14) 


(2. 15) 


This transformation allows f and to be unbounded. 

In order to interface the high and low thrust stages a correspondence or 
scale factor is needed to relate the respective costates where 


htf 


c = 




HIGH 


and t^ and refer to the respective costs for the two stages. Then 




LOW 


HIGH 


(2. 16) 


(2. 17) 


The iteration parameters are then 


X - 


^2 

^3 


c 

fi or to 
Xn 


(2. 18) 
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and recall that the initial a, e, i, and o) oi" Q are specified. 

If roll is free so that the panels can always face the sun and if shadowing 
is not included, then the perturbations on the spacecraft are axially symmetric. 
The perturbations modeled are oblateness and the Van Allen radiation. The latter 
was forced to be axial symmetric by the particular model used. Therefore, X 
is theoretically constant and varying O at the initial time will only affect q at 
the final time. Thus if p at the final time is not specified, the sensitivity matrix 
will be singular. For this case the q. should be specified. 

The sensitivity matrix is calculated by varying each x. by a small amount 
and comparing the resulting error in the final conditions with the nominal. 


n 




(2. 19) 


ax. 


6X; 


Except for — the partial of y with respect to t^ can be calculated analytically 


3 tf 

since the derivatives of the state and costate are known at the final time, 
example 


For 


dYo 

— = h(t ) (2.20) 

The other partials are obtained similarly except for the partial of the Hamiltonian. 

For SEP the time dependence of the Hamiltonian derives from the shadow 
effect and solar distance and direction and therefore an analytical evaluation 
is difficult. Therefore, the partial was obtained numerically by varying t^ 
slightly and re-evaluating the state derivative (holding the state and costate constant) 
and then evaluating the Hamiltonian. Thus 


^H(t^) X^^(t^) ^x(t^+^ t^) - x(t^)^t^^ 

^ 

^t|. ^t^ 


( 2 . 21 ) 


2.4 Numerical Methods 

A Newton-Raphson iterator is used which calculates the sensitivity matrix 
by running neighboring trajectories by changing slightly the initial values of the 
iteration parameters, one at a time. The size of the change in the iteration 
variables is chosen by the user and can affect the accuracy of the matrix. A 
Modified Newton-Raphson iterator uses basically the same technique but many of 
the iterations make use of a modified sensitivity matrix rather than calculating a 
new one by running neighboring trajectories at each iteration. 
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The low thrust differential equations are integrated using a fourth order 
Hunge-Kutta method. The time step is selected by the user. This is fixed for 
one subroutine supplied or may internally be changed if the IBM Scientific Subroutine 
version is in which case error weights and an upper error bound must be 

supplied by the user. Cutting the size of the time step can increase the accuracy 
of the trajectory but rapidly increase run time. 

Numerical averaging utilizes a Gaussian quadrature. The number of points 
sampled on an orbit can largely be determined by the user. Again^more points 
increase accuracy at the expense of run time. 
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SECTION 3 


LOW THRUST EQUATIONS 


3. 1 Introduction 

This section developes most of the new analytic results. Some of these new 
results are based on previous work, which is included in summary form, often 
in appendices. The new results include the derivation of the attitude constrained 
thrust control, the effect of a delay in thruster turn-on after leaving the shadow, 
the new solar array degradation model and the calculation of several spacecraft 
parameters. 

3.2 Equinoctial Orbital Elements 

By using equinoctial orbital elements the singularities that occur for zero 
eccentricity or inclinations of zero or ninety degrees when using classical orbital 
elements are avoided. (For inclinations near 180^ retrograde equinoctial orbital 
elements can be used, although we will not consider that case in this report. ) The 
formulas given in this section are taken from Reference 16 except for the costate 
transformations. 

The direct equinoctial orbital elements are defined in terms of the classical 
orbital elements, a, e, i, p , and by the formulas 

a = a 

h = esin((o + 0) 

k = e cos (o) + 0 ) (3.1) 

p = tan (.j) sin o 
q = tan (-j) cos Q 

In Reference 16 the sixth orbital element is the mean longitude at epoch. In this 
paper we will consider the eccentric longitude, F, as the sixth element, defined by 

F = E + CO + p (3.2) 
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where E Is the eccentric anomaly. This element will be eliminated from the 
dynamical equations by the averaging process. 

The inverse relationships are defined by 


a 

e 

i 

Q 

u) 


tan ^(p/q) 

tan ^(h/k) - tan ^(p/q) 


(3.3) 


The equinoctial coordinate frame is defined by the basis vectors g , w , 
which are given below with respect to an earth equatorial coordinate frame. 


f = 


5 2 " 

1 + p +q 


1 2 ^ 2 - 
1 - p + q 

2pq 

-2pq 


I 5 2 " 

1 + P +q 


w - 


1+7 + q 


2pq 

1^2 2 
1 + p - q 

L 2q 

2p 
-2q 

II 2 2 

|_1 - p - q . 


(3.4) 


This coordinate frame is illustrated in Figure 3-1 where w is normal to the 
orbital plane. 

The equinoctial orbital elements can be calculated from position and 
velocity. The semimajor axis is given by 


a = 



(3.5) 
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From the values of p and q and Eq. (3, 4) £ and g may be computed. Then 

h = e • g 
k = e . f 


(3. 10) 


Further relationships include the position coordinates Xj, Yj with respect 
to the f, g, w frame. 


Xj = £ • L 


Yj = £ • I 


(3. 11) 


and the eccentric longitude, F, where 


cos F = k + 


(1 - k^,5)Xj - hk^Yj 


ijl - h^ - k^ 


(3. 12) 


sin F = h + 


(1 - k'=0)Yj - hk^Xj 


1,/! - h - k 


where 


B 


(3. 13) 


i./TTiTT?- 

The mean longitude is defined by 

X = M + o) + n 


(3. 14) 


The eccentric longitude, F, was defined in Eq. (3.2). Kepler's equation in terms 
of X and F is then given by 


X = F - ksinF + hcosF 
We will not make use of the true longitude 

L = V + to + O 


(3. 15) 


(3. 16) 
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M, E, V are the mean, eccentric, and true anomalies, respectively. 
Position and velocity are given by 

r = Xjf+Yjg 

, * ^ * A 

r = 

where 

= a[(l-h^^)cos F + hk|9sinF - k] 

Yj = a[{l-k^ g ) sinF hk/? cos F - h] 

- ng.2 2 

[hkgcos F - (1-h fl) sin F] 

• na^ 2 

[ (1-k^^ ) cos F - hk^ sin F] 

and 

— = 1 - kcosF - hsLnF 
a 


( 3 . 17 ) 
( 3 . 18 ) 

( 3 . 19 ) 
( 3 , 20 ) 

( 3 . 21 ) 

( 3 . 22 ) 

( 3 . 23 ) 

( 3 . 24 ) 


fi is the earth gravitational constant. 

The adjoints to the classical and the equinoctial orbital elements are related. 
Let X represent the adjoints to a, h, k, p and q, and let ^ represent the adjoints to 


a, e, i, n and Then 


'^a' 


" 1 0 

^h 


0 sin(co+n) 

>^k 

= 

0 cos(to+n) 

Xp 


0 0 

Xq 


0 0 


0 

0 

0 

0 

0 

0 

2 i 

2sinQcos 

cosp / tan -5^ 

o 2i 

2cosp cos ^ 

-Sinn /tan^ 


0 


t 

I 

cos(tri+0)/e 



-sin((o +0 )/e 



-COSO /tan^ 



sinp /tan-g 


60 


( 3 . 25 ) 
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This can be written Ln terms of equinoctial orbital elements 



Note that some of the elements of this matrix are undefined if e or i are zero. 

The inverse relationship can also be written. In terms of classical 
orbital elements, the classical costate is given in terms of the equinoctial costate 
by the following relation. 



1 0 

0 sin(^+Q) 

0 0 

0 ecos(to+o) 

0 ecos(t^+Q) 


0 

cos(u3+0 ) 

0 

-esin((*)+Q ) 
-esin(o;^+Q ) 


0 

0 

1 2 i . 

•jsec ^sinn 

tan-jcosp 

0 



( 3 . 27 ) 


or in terms of equinoctial orbital elements 



( 3 . 28 ) 
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Certain elements are undefined if e or i are zero. 


3. 3 The State Equations 

The unaveraged state differential equations are 

2 = — M(z, F) 

““ me 


m 


2P 

c 


(3. 29) 
(3. 30) 


where z represents the five equinoctial orbital elements (a, h, k, p, q), m is 
spacecraft mass, P is thruster beam power, which is assumed given by 


P 


T]P D(N) 
o 



cos 0 


(3.31) 


where P^ is a constant initial maximum array power at 1 A.U. , rf is the total 
constant power efficiency factor, D(N) is a damage factor which is a function of the 
fluence, N, and has a maximum value of 1 (this function is derived in Section 4*3), 
is the distance to the sun in A.U. and 0 is the angle that the normal to the 
panels makes with the spacecraft - sun vector. For the unconstrained attitude case 
cos 0 = 1. For the attitude constrained case it is a function of spacecraft position 
and the thrust direction. The constant exhaust velocity is given by 


c = 


sp 


(3.32) 


where I is the specific impulse and g^ is the acceleration of gravity at the equator, 
sp 0 

The thrust direction is given by the unit vector u . We define the 5x3 matrix 


Bz 

M(z,F) = — (3.33) 

The elements of this matrix are listed in Table B-1 of Appendix B. In addition to 
the quantities defined in Section 3. 2 the partials of and with respect to h and 
k are required. These are listed in Table B-2. These partials differ from 
Reference 16 since we consider F as an orbital element rather than as a function of 
h and k; thus when partials are taken, F is held constant. This assumption also 
affects the appearance of the expressions given in Table B-1 when compared with 
Table 3 of Reference 16. However, if the expressions are written out in detail 
they are seen to be identical. The code utilized the form shown in Ref, 16. 
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That form of M was also used when calculating and coding ^ — for the costate 

(2 8 ) 

equations. An algebraically simpler form of this matrix is given in a later work' 
but was not utilized in the coding of our program, although a variation of it is 
discussed in Appendix C. 

The high energy particle fluence is also included as the seventh state variable. 
It is a function of position, but the exact form is developed in Section 4.2. For 
the present discussion assume 

N = f(z, F) (3.34) 

When mass flow rate is given by Eq. (3.30) and exhaust velocity is constant, 
for arbitrary power history, the aV is given by 

m . 

-c In (— 5 -) (3.35) 

V(t)^ 


3.4 The Canonical Equations 

To simplify the analysis we will first consider only the costate elements 

z, m. The effect of fluence can be considered later. For simplicity also assume 

that in Eq. (3.31) that t) = 1 and R = 1. Then beam power is given by 

' s 

P = P cos 0 (3,36) 

o 

The method used to get the canonical equations is to form a Hamiltonian, 

H, then average it, from which a first order approximation to the actual state and 
costate is obtained. (This does not necessarily yield the averaged costate. ) 

Let H represent the unaveraged Hamiltonian and x the full state, then the 
averaged Hamiltonian is 

T 

t+-f ^ 

= — ^ T h(x X.r.FfU.ufFfOl^dt = ^ ^ H(x,X.t.F.u(F))-^dF 

T ^ ^ T^-tt^ dF 

o t- -2“ o 

(3. 37) 

where T^ is the orbital period, and the overbar indicates approximate quantities. 
The F dependence refers to time dependence which is independent of spacecraft 
motion such as Earth motion around the sun. In practice integration with respect 
to F is more convenient. The expression can be obtained from the equinoctial 
version of Kepler*s equation (Eq. (3. 15)). It is given by 
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T 

Hi = _£ (1 - IfcosF - hsinF) (3,38) 

dF 2ir 

For convenience define 

s(z.F) = — ^ (3.39) 

T dF 
o 

Then the canonical equations for the approximation to the state and costate are given 

by 


I = ^ = r ^ s dF 
aX 

i = = -r \^s + hmI)cif 

ax ax sx 


(3.40) 


(3.41) 


Note that inside the averaging integral the state and costate are held constant. From 
the definition of s the only non-zero elements of ^ are 

c^x 


SS 1 . T-. 

— = - — sinF 

2tr 

^ = - -1 cos F (3.43) 

ak 2 it 

is just the unaveraged state equations with the dependence on the state and 
B X -V H 

costate held constant over the averaging interval. The expression is just 

• ^ 2L 

the unaveraged costate equation, x,- Most of the remainder of this section will be 
involved in obtaining the unaveraged Hamiltonian and its partial derivative with 
respect to the state for the attitude constrained case. Those equations reduce 
simply for the pitch constrained or unconstrained cases. Only the state variables 
z and m will be included; shadowing and oblateness will not be included. 

One method of including the constraints is to write u directly as a function 
of a single control angle, which lies in a plane perpendicular to the radius 
vector. Let 


T 


1 



(3.44) 


22 



where and 62 are any two unit vectors orthogonal to the radius vector and each 
other (the transformation Is thus dependent on the spacecraft position and thus 
the orbital elements). Then 


z = 2P{z,V.^) M(z,F ) Tj(z,F)[^°®^] 


(3. 45) 


rii 


“ — 2” F , i/} ) 

c 


(3. 46) 


The Hamiltonian is given by 


H 


.T 2P 
— z 

me 



- X 


m 


2P 

2 


The costate equations are 


X 

— z 


li - — I— T 

P -z me ' 1 


+ M — - \ 

az ' 



(3.47) 


(3. 48) 




—2 


2P 
m c 



(3. 49) 


The elements of expression are given in Tables B -3 to B- 7 of Appendix B along 

^z 

with subsidiary partials in Tables B -8 to B-9. 

Another method is to use a constraint formulation (which is the form coded). 


Then 


2P 


M u 


me — 


m 


2P 

c 


but there is a constraint equation 


C 


T 


r u = 


0 


(3. 50) 
(3.51) 


(3. 52) 


which insures that the thrust direction is orthogonal to the spacecraft radius vector. 
The Hamiltonian is 


H 


.T 2P ^ „ T- 2P . 

X Mu + ^r u - — jc X 

-z me - - 2 m 

c 


(3. 53) 
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with the additional multiplier, U . Then the costate equations are 


32 -z 


2 

Mu ^ 

me - 

+ X 

2P 

(BM + y 911 1 , Q 

(3.54) 

J “Z 

me 

' gz j - 

II 

Mu 



(3. 55) 


AeaLn. Is given in Appendix B. If we assume that r is the radius vector divided 

- SXj SYi SYi 

by the semimajor axis then the only nonzero parUals of r are ^ ’ 

which are given in Table B-2. Thus the term is a little easier to evaluate than 
- 2 — , and also this form is more useful for the cases of zero pitch but free roll and 
\vithout constraints. However, ll must be evaluated. Rewrite the Hamiltonian 


(\J ^ M + fir ■ 
me 


^ ^ j- \ 

U “ — 77 A 

z m 


= (xj — M + fi r'^) Tt"^ £ ■ \ 


(3. 56) 


where 




(3.57) 


(3. 58) 


Now T^^u is just the control in the , £2 * £3 coordinate system. It can be written 


cos $ cos ^ 
sin il) cos 
sin ^ 


(3. 59) 


where p is the angle between u,p and the e^ - £2 plane, and is the angle between 
the projection of u^ onto the “ £2 have not 

taken into account the constraint equation. Rewriting Eq. 3.56 


me c 


(3. 60) 


Writing this out in detail 


op Tf* A A P2P T A n 2P 

H = X M[e j cos 0 + £2 sin 4) ] cos ^ + I — X^ Me^ - »i |r j j sin f - 

me me e 


(3.61) 
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The longitudinal axis of the panels Is parallel to the pitch axis. Therefore, 
the orientation of the panels can be defined independently of the pitch angle. Thus 
P is not a function of pitch. Then to solve for the optimal f set 

ML = 0 (3.G2) 


or 


tan ^ 


,/ 2P . T 
Vmc -z 


M 


-3 



^ M[e j cos If) + £2 sin ^ ] 


(3. 63) 


Now the constraint can be taken into account. Rewriting Eq. (3. 52) as 


C = r”^ Tt"^ u = 0 (3. 64) 

since Tt"^ is just the identity matrix. From Eqs. (3.57) and (3.59) the constraint 
reduces to sin f = 0. So Eq. (3. 63) is also zero. Thus 


2P 

me 


Me„ 
Z — O 


Mir I 


0 


(3. 65) 


or (using Eq. (3.58)) 


4 


2P 

me 





(3. 66) 


When roll is free, power is no longer a function of $. For the pitch = 0 
case Eq. (3.52) is still valid, and the optimal u is easy to calculate as is . The 
Hamiltonian is 


H = (XT^P M + 4£^)u-^ X^ 


(3. 67) 


me 


Since P is not a function of fi , to maximize H set 


u - 


(X 


z me / 


ix"^ — M + /ir'^1 
1 -^ me > 


(3.68) 
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The constraint C = u r = 0 yields 


T 2P T 

X — M r + Mr r 
— z me — — - 


= 0 


or 


me Irl'^ 


the same as Eq. (3. 66). 

When there Ls no constraint n = 0 and 


T 




is the optimal thrust direction. 

In summary the approximating canonical equations are 


z = ^ f2 P(|, u) vi(z,F) + 4r’^)u s(z.F)dF 

7T ^ E^iC 

. p IT 2P(z,u) _ 

m =-\ 2 s(z, F)dF 

'^-1T c 


i . - [ "> Mg Fin -■L)*T 2P(i,uY8M(£.F ) 

^ 3z ^ me c ^ me , Bz 

+ fi ^-»^^ ^u|s(z,F) + ]dF 


* i r* ^ 

X = — \ >• z(z,X.u)s(z.F)dF 

m m J_„ -z 


where 


^ . _jT 2p(;,j) „g )_£ 

— z — — 


me 


l£l 
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(3. 70) 


(3. 81) 


(3. 72) 
(3. 73) 


(3. 74) 


(3. 75) 


(3. 76) 
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if pitch is constrained to zero. If there are no constraints, tt = 0. M and — 

az 

are given in Appendix B. Also 


0 


as ^ 1 

az 2 tt 


- sin F 

- cos F 

0 

0 


(3. 77) 


If r is Interpreted as the radius vector divided by the semimajor axis, then in 
the equinoctial coordinate frame 


so that 


r 


X^/aT 

Yj/a 



0 0 0 1 


^1/a 


ah " 


0 


^^1/a 

at 


^Yl/a 

ak 


0 


0 0 0 
0 0 0 


(3. 78) 


(3. 79) 


where the nonzero partials were given in Table 

aP 

In the case of free roll, P is not a function of z and so — =0. For the 
zero pitch and roll case this partial depends on the actual form of P which has 
not yet been specified. For the unconstrained case and for the case with zero pitch 
but free roll, the canonical equations are completely defined. It remains to 
calculate the power equation and the control for the zero pitch and roll case. 


3 , 5 Geometry and Power 

In this section we define a useful coordinate system in which it is convenient 
to calculate the thrust vector direction and the equation for power. It is shown that 
we can calculate the panel orientation angle in terms of the thrust direction angle. 
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thus leaving only one control variable for the problem. Power will then be a function 
of the thrust angle and the angle between the radius vector and the direction to the 
sun. 

The thrust vector Is constrained to a plane which Is perpendicular to the 
radius vector. In this plane the yaw is usually measured from a basis vector in the 
orbital plane, which it is for the printed output of our program. However, because 
of the power dependence on sun direction, another coordinate system will be useful. 
This is the e system where 

e 3 = - r (3.80) 

and e ^ is in the plane defined by the radius vector and the Earth-sun vector 
(assumed equivalent to the spacecraft-sun vector). In particular. 


and 


e^ - Og X £2 

In this system the sun direction is defined by the angle ^ where 

= cos 63 + sin p e j (3.83) 

^ T A 

(In the equinoctial system cos^ ' " £ 5s * ^ thrust acceleration vector direction 

is constrained to the e ^ - e^ plane and given in terms of the angle ip with respect to 
the e^^ axis by the unit vector 

u “ cos + sin (3.84) 


(See Fig. 3-2). 

The power yield of the solar panels is proportional to the cosine of the angle 
between the normal to the panels and the vector pointing toward the sun. The panels 
can rotate about an axis which is in the e^ - £2 plane and is assumed to be perpen- 
dicular to u (recall Fig. 1-1). The normal to the panels can be defined in terms 
of ij) and an angle of rotation, 0, about this axis in the " i-2 “ ~3 namely: 

N = sin0cos0£^ + sin0sin0£2 + 003063 (3.85) 
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FLg. 3-2 Coordinate and Angle Definition 

The cosine of the angle between the normal and the sun vector is given by 

COS0 = = sLnjSsLne cos»(^ + cosecos^ (3. 86) 


For a given q it is possible to maximize this expression with respect to 0 and so 
eliminate 0 from cos^, which will then be a function only of ij) and The 
maximizing 0 is given by 


sin0 = 


sin S cos ^ 

/ 2 !~2 2 

ycos -+* sin flcos 0 


COS0 = 


cos B 


y o 2 ' 

cos^g + sin gcos^/> 


(3. 87) 


(3, 88) 


and thus 


COS0 




"2 2 2 
COS g + sin geos ip 


-fi 


2 — 7 ~ 2 — 

sin psin ij) 


(3. 89) 
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It is interesting to look in the e^ - plane at the locus of the vector pointing 
along the thrust direction with magnitude equal to cos 0. For each £? there is a curve 
in the plane as varies from 0 to 360 degrees. These curves are symmetric 

about both axes and repeat for every p interval of 90 degrees. Note that for |cos/? ] 

< yT/2 they are convex (see Fig. 3-3). Since power is a function of thrust direction, 
this curve will be important in determining the desired control angle, which 
specifies the thrust direction. 



Fig. 3-3 P/Pq Locus 


In Eq. (3. 74) for the costate derivative we needed . This can now be 

oZ 


calculated. Power is given by 


P = P cos 0 = P Jl - sin psin fp 


(3.90) 
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so that 


^ = - P (l - sin^/^sin^i/)) ^ sm^l()SLn^cosp 

o\ ^ ) o£ 


(3.91) 


But cos^ ~ • Rg . So 


where 


iSSli . . R + e, . ^ 

az 9Z s -3 az 


(3, 92) 


3£3 

aZ; 


III 


V’'? Pf - Vi -R-) 

/ 2 \l 


(3. 93) 


where the semimajor axis drops out of and so we only need the partials of 
and Yj with respect to h and k which are given in Table B-2 of Appendix B. 
is a function of only p and q and the partials are given in Appendix D. 


3. 6 Thrust Direction 

Let us rewrite the Hamiltonian given in Eq, (3.47) but including dependence 
of P on i/) . 


H = 


2P 

o 


I ! 5 2 

^1-sin psin ij) 


T MT 
m 




(3. 94) 


In the definition of the transformation, T, we are using the unit vectors defined 
previously in Eqs. (3. 80). (3. 81), (3. 82). The quantity is just the primer 

vector given in the system and will be designated ^ with components Xv|*^v2' 

Xv 3 « Eq. (3. 94) can be rewritten as 


H = 


yi-sln^gsln^ij) (^-^cosi/) +-g^sini^) - ^) (3.95) 
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For a minimum time solution this Hamiltonian must be maxLmizeci with respect 
to 0. If this is done by setting = 0 and converting all sin 4) to cos 4> , a 

sixth order polynomial in cos ip results. Let 


e 


mX 


m 



tan^ 



(3. 97) 


Then 


H 


2P 


cm 


X ^/cos'^jfj + sin^^cos^^ (cos0 + tan^ysin^ - c) 


(3.98) 


This can be maximized with respect to ip by setting 


^ = 0 
dll) 


(3. 99) 


u 2P 
aH ^ o 

dib cm 


■ 1/2 


r/2 2 2\ ^ 2 / N 

X^ -{“^cos ^sLn jQcos sin |5cos4^ sini/; ^cos0+tanQsinj|>+€ j 


^l/2 


+ ^cos^l9+sin^j3cos^j^^ ^-sin^/>+tan(y cos0^ = 0 


(3. 100) 


2P X 

oV, 


Divide by 


2 2 2 1 /2 
cm(cos jP+sin pcos ’ 


, then 


or 


2 2 2 2 
sin pcos0sint() (cosif)+tanasinj/j+f ) + (cos g+sin f?cos if) K-sin^j-ftanc^cosj/) ) = 0 

(3. 101) 


2 2 2 2 2 2 
“ sin pcos if)sinif) - sin ptainotcosipsin ip - csin jgcosi/)sini/) - cos gsinip 

2 2 2 2 3 

+ cos jgtanacos^ - sin ^gcos i/>sinif) + sin i^tan^cos i/) = 0 (3.102) 


Rearranging 


2222 22 23 

sLnif)(2sLn geos 0+cos ^ + csin gcosif)) = (cos j^-sin g)tanacosi|)+2sin ^tan^cos ip 

(3. 103) 
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Squaring 


(1-cos 

. 2 . 

+ e sin 


2 

4 


0 )(4sin'^gcos^if)+4sin^pcos^ | 9 cos^iIj+cos 0+2esin ^(2sin jjcos 0+cos p) 

2 2 2222 22 2 2 4 

= (cos p-sin p) tan c^cos 0+4sin 0tan ^(cos g-sin p)cos 0 

+ 4sin^^tan^fvcos^0 (3. 104) 


Rearranging, finally. 


4 2 

4sin flftan r 


.4-1 \ r\ c? 


+ 4rsin^ p tan^o (cos^ 0 -sin^ 0)-stn‘^ft+sin^0 cos^fl - e^sin^g] cos^i/j 

4 2 2 3 

+ e(4sin |9 - 2sin pcos p)cos 0 


r 2 222 2 2 4 2.4n2 

+ (cos | 5 -sin tan ev^sin pcos ^+cos p+c sin ft | cos it 

2 2 4 

+ 2c sin pcos Qcos^ - cos p = 0 


(3. 105) 


Is the minimum time solution really desirable for this problem? We think 

not. If power were not a function of the control, the optimal thrust direction 

would just line up with the projection of the primer vector on the e^ - £2 Plane. 

But for the Hamiltonian of Eq. (3. 95) 0 becomes a function of This has the 

effect of biasing the control so that the thruster is operating in a region on the 

curves of Fig. 3-3 where cos0 is greater than it would be were not included. 

This reduces the mass of the spacecraft by throwing fuel away in order to reduce 

the flight time. In fact, if fuel could simply be dumped overboard instead of going 

through the thruster, the time optimal requirement would be to dump fuel until 

X =0. Since the formulation here does not allow separate dumping of mass, a 
m 

bias on the control is required to reduce mass. Although minimizing time is 
important, dumping fuel to do so seems undesirable. 

So what happens if the part of the Hamiltonian containing X^ is ignored and 
instead we maximize 


H' 



(X COS0 + X sin0) 

'^2 


■? 


(3. 106) 


In the 6j - eg plane one can see that this is equivalent to maximizing the projection 
of the vector (cosip, sini/; >A -sln^g sin'll/) onto the e^ - ig components of the primer 
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, Xv 2 )* locus of the first vector Is just that shown Ln Fig. 3-3. Geomet- 

rically, a typical solution is shown In Fig. 3-4. The maximization of H' can be 
carried out using analytic geometry by considering the cos^ curve as a parametric 
function of ij) and noting that its slope must be the negative of the inverse of the 
slope of the Xy^ - X^^ vector. Thus, If 


x(!/j ) 


cosli) 

y(il)) 


sin0 


yi-sin^flsin 


(3. 107) 


then 


dy X 

^ ^ = - 1 = ^ 

dx to tanft X 

dll) '^2 


(3. 108) 



Fig. 3-4 Geometric Calculation of j/i 


where 


2 2 2 2 

to „ (-1+sln jSsin ij)-sin geos l/))sini/) 

dip 2 2 ^ 

(1-sin psln^l^) 


(3. 109) 


34 


and 


2 2 

dy „ (l-2sLn ffsin (/j)cos^ 

dip 2 ~ 2 ^ 

(1-sin flsin^0) 

Thus 


Using 


and 


then 


2 2 2 2 

s _ 1-sin 0sin 0+sin flcos ib . , 

tau^ - ^ ^ — 2 tan;/) 

l-2sin flsin ip 


. 2 
sin 


0 


. 2 , 
tan w 

1+tan 


2 , 

cos tp 


T+tan^^ 


(3. 110) 


(3. Ill) 


(3. 112) 


(3. 113) 


tanft 


2 2 2 
l+sin ft+cos jgtan 0 

1+tanjf) -2sin 5 tan ip 


tan$ 


(3. 114) 


or 


cos^0tan^0 + (2sin^f?-l)tana?tan^4> + (1+sin^ft )tani/) - tan« = 0 (3.115) 

2 

Dividing through by cos ft yields 

tan^i^ + (tan^/9-l)tan«tan^i/j + (2tan^ “f-l)tanj/) - (1+tan^fi )tanev = 0 (3.116) 

It is also possible to obtain the control by setting = 0, The procesure 

d0 

is the same as resulted in Eq. (3. 105) and in fact an equation in cosip can be 
obtained directly from Eq. (3. 105) by setting e ^ 0. Thus the coefficients of all 
odd powers of cost/) are zero. The resulting equation is 

4sin^j3(tan^a + l)cos^(/j+4sin^ P(tan^«“ IKcos^^ - sin^ ^?)cos^0 

r/2 222 2 2 4 2 4 

+1 (cos )9-sin p) tan <j-4sin ficos p+cos ^)cos ip-cos = 0 (3.117) 
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The cubic equation in tan^ of Eq. (3.115) can be shown to be equivalent 


to the cubic in of Eq. (3. 117) by substituting tanj/) = - — into 

cos 

Eq. (3. 115), then removing the square root, rearranging and squaring. Collecting 
the terms and making certain trigonmetric identifications yields Eq. (3.117). 

This control does not dump fuel to minimize time, and it is much easier to 
calculate since the solution to the cubic can be obtained analytically, unlike the 
solution to the sixth order polynomial. Analysis of the worst case example (see 
Appendix E) indicates that the time penalty in using this suboptimal control is 
small while resulting in a fuel savings. Since it is also much easier to calculate 
and thus saves computer time, it was decided to use this control law for the present 
study. 

This control law has some interesting characteristics. Fig. 3~5 shows 
curves of constant « superimposed on the curves of Fig. 3-3 for one quadrant. 

Thus for a given p (sun angle) and ct (primer vector angle) the resulting is 
just the angle made by the line from the origin to the intersection of the appropriate 
ft and p curves. Note that in this quadrant ip is always less than ft except for 
0=0 when ip = a or for p < 45® and ft = 90® when also ip = 90®; when 0 = 90® , 
ip — 45® as ft 90® . For 0 > 45® as ft crosses the e 2 axis there is a jump in ip 
to a symmetrical position in the adjoining quadrant. If o remained aligned with 
the £2 axis for a finite time a singular arc results. Since it would be nonoptimizing 
(in terms of HM to operate on the concave portions of the curves, a chattering 
solution is required as i() jumps back and forth in the two quadrants in infinitesimal 
time yielding a resultant thrust vector along the £2 axis. In practice this solution 
will probably never be required (as well as being impossible to implement if it 
were). Jumps in the control direction can and do occur however. The locus of 
points on Fig. 3-5 at which a jump occurs is given by an arc of a circle of radius 
J^l2 and center at the origin. This radius represents the minimum power. The 
at which a jump occurs is related to by the relation 

COS^^ip = ~n ^ (3.118) 

sin^0 

In practice, then, 0 and ft are calculated and the coefficients of Eq. (3. 117) 
are obtained. The cubic equation can then be solved analytically. Appropriate 
modifications must be made if the tangents are very large or infinite. 
There are three roots to the cubic equation and if all three are real 
they correspond to six possible values of in a 360® range. These 
correspond to various maxima and minima. By inspection of the form of 
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the figures produced by Eq. (3. 107) there is apparently only one solution which 
can lie in the same quadrant as and at the same time maximize the projection 
onto (Xvi * Xv 2 ^- Thus we can obtain the control angle tp which can then be 
substituted into the equations of motion. This calculation must be made at each 
quadrature point on an orbit. 





Fig. 3-5 Lines of Constant Primer Angle (o? ) 
and Sun Angle (g) 

There are certain special cases of values of a and g that must be considered. 
If tanaf = 0 the cubic reduces to tan^ii) + (2tan^ « + l)tan?i) = 0. The correct root is 
then tan^ = 0 and thus ili = a- If tana is infinite, first divide Eq. (3. 115) through 
by tana , then set tan« equal to infinity, which yields 

(2sin^0-l)tan^ii) -1 = 0 (3. 119) 
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or 


or 


tan 0 = 


2sin 


2 ; 

COS 0 


2sin^fl-l 
Sin R 


(3. 12 0) 


(3. 12 1) 


2 1 

This is valid if sin /? > ^ . Otherwise tan0 equal to infinity is the correct root 
of Eq. (3. 115) and ip = a (- 90° or 270^). If cos^^ = 0, then Eq. (3. 115) reduces to 


2 

tanfttan ip + 2tan^/) - tana = 0 
From the binomial theorem 


(3. 122) 


tan0 - — ^ ±yi"tan a ^ tan i o , - cot ~ « (3.123 

tanfv ^ 

If -90° < a < 90° , the first root is correct and ^ = i if 90° < a < 270° the 

1 " 1 * 
second root is correct and i/j = -f 90° if 90° < a < 180° and ih = - ot - 90° if 

2 ^ ^2 
180° < a < 270°. If cos 0-1 then tani/) = tana is the correct root and = a. 

Because of the numerical inaccuracy of using a quadrature over a region 

which contains a discontinuity in the integrand, it is best to calculate possible 

jump points and then do separate quadratures be tween these points. A jump in the 

thrust direction occurs when the projection of the primer vector onto the 

^2 ” —2 passes from one side of lo the other. At the time of crossing 

the projection is perpendicular to e , and also perpendicular to R . Thus the 
• — ® 
condition for a jump can be written 


. T T T 



= 0 


(3. 124) 


An actual jump occurs only if the constant curve is concave, i.e., |cosj9 [ , 

For the results given in this report the method used to find the values of F for 
which a jump might occur was to divide the 360° range of F into small equal 
intervals and then to check to see if there was a sign change of Eq. (3. 124) from 
one side of an interval to the other side. If so, the exact F within the interval 
for which Eq. (3. 124/ was zero was found by using a nonlinear function root 
finding routine (an IBM scientific subroutine using a Mueller's iteration^^"^\ 
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The equation for the jump points can be obtained as a sixth order poly- 
nomial in cosF. This equation is very tedious to derive algebraically and in any 
case must be solved by an iterative procedure and the roots tested. Thus it seems 
reasonable to use the search procedure. The sixth order polynomial in cosF can 
be obtained by writing a form of Eq. (3. 124) in r, s, w coordinates where r is 
along the spacecraft radius vector, w is perpendicular to the orbit plane and 
s = w X r . If X and X are the components In the s and _w directions respectively, 

then 




+ «v-rJ>^s -z-w 


+ R x„ 


= 0 


(3. 125) 


where = (R , R * R^) the sun's direction in equinoctial coordinates and 
and are the spacecraft position in the equinoctial frame, and r is the radial 
distance. Using the results of Appendix C 


= 2a7l-h^-k^ 
^s ' nr 

G 




2 , 2 > 


(qY^-pX,) t (XpY,n^X,) 


(3. 126) 

(3. 127) 


The quantities X^. Y^, X^. Y^. r all contain cos F and sin F. Substituting into 
Eq. (3. 125), multiplying through by r and eliminating sinF by rearranging and 
squaring results in the sixth order polynomial in cosF. Three of the roots are 
introduced by squaring so there will be at most 6 values of F for which there may 
be jumps. 

If the eccentricity is zero, h = k = 0, and a quartic results. In that case 


Xg = 2 (aX^ + Xj^sinF + XqCOsF\ (3. 128) 

= /T (X SinF + X^cosF) (3. 129) 

W ^ ^ ^ P H 

There are at most four values of F at which there may be jumps. If in addition 

to h = k = 0 also X, = X, = 0, then there is a further simplification and Eq. (3. 125) 
Q K 

reduces to 


tanF 




(3. 130) 
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where 


R (l+p^+q^) 

K = — 

4aXa 


(3. 131) 


Thus there can be at most two jumps which are 180° apart for this case. 


3.7 Qblateness 

In the previous sections we have considered only perturbations to the 
inverse square motion caused by thrusting. In this section the effect of oblateness 
(J 2 ) is considered. Oblateness was considered in Kef. 2 and the equations are 
included here for convenience. Oblateness does not directly enter the power and 
mass derivatives so only its effect on the orbital element derivatives will be shown 
in this section. 


The single-averaged perturbing potential due to Tg has been calculated in 
terms of equinoctial coordinates in Ref. 16 and is repeated in Appendix F. R^ is 
the equatorial radius of the earth and J 2 is set to .001827. These formulas enter 
the averaged Hamiltonian as coefficients of the costate (outside the integral since 
the averaging effect has already been accounted for). 

If ~z indicates the perturbation due to thrust as given in Eq. (3. 72), 

ci 

then the Hamiltonian is given by 

H = x"^z - + X% (3. 132) 

J2 - -a 

The state equation Is 

z = z , + z (3. 133) 

^ “”0 rt ' "di 


The costate equation is 




az az -fr *- gz ® SZ 


(3. 134) 


The partials indicated by , /^z in the above expression are given in Appendix F. 

^2 ” 
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3.8 The Shadow Effect 


For SEP missions, the thrusting will be shut off while the spacecraft is in 
the earth's shadow. The entry and exit angles are needed in order to perform the 
averaging integral. In calculating these angles the following assumptions are made. 
The shadow is cylindrical; the earth revolves around the sun in an elliptical orbit; 
and over one spacecraft revolution, the sun's direction is fixed. The time of 
thruster turn-on can either be immediately upon exit from the shadow or after a 
delay following the exit from the shadow. The delay time calculation will be 
considered in the next section. Calculation of the entry and exit eccentric longitudes 
was considered in Ref. 2. The pertinent equations are summarized in Appendix G. 

Let F refer to the eccentric longitude at entry to the shadow and F^ at 
exit (or the thruster turn-on time if a delay is included). That part of the Hamiltonian 
proportional to thrust is then 

H = 


Thus 


X - 


H s dF 
F, 


2 

X s dF 


(3. 135) 


(3. 136) 


and by Leibnitz' rule 



(3, 137) 


The calculation of ^ is discussed in Appendix G and the modification needed if 
Includes a delay"effect Ls discussed in the next section. 


3.9 Delay in Thruster Start-Up After Leaving Shadow 

Some orbits of the spacecraft may pass through the shadow of the earth. 
During this period the thrusters are assumed to be turned off. The shadow entry 
and exit points are derived as the solution of a quartic equation and given in 
Appendix G under the assumptions of a cylindrical shadow and stationary sun over 
one orbital period. It was assumed in Ref. 2 that the thrusters were turned on 
immediately upon leaving the shadow. A more accurate model is to include a delay 
in turn on time. This delay time is the sum of the time for the solar array to 
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achieve operating temperature and the time for the thruster to achieve full thrust 
after the solar array power is applied to the power processor. 

The model used is taken from Ref. 2 9, The thruster start-up time is 
modeled as a quadratic function of time in shadow. The solar array temperature 
adjustment is modeled as a constant of 2 minutes. The total delay is then just 


= aT ^ + br + c 

where is the time in shadow. If t and are in minutes, then 

1 

^ Y7UU 

' 77UU 

c = 12 


(3. 138) 


(3. 139) 


This function is plotted in BTg. 3-6. 



Fig. 3-6 Delay in Thruster Turn-On as Function 
of Time in Shadow 
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The shadow entry and exit eccentric longitudes are known. From these the 
time In shadow is obtained from Kepler's equation. In equinoctial elements the 
time In shadow is 

^ ^ [(Fex - Fen) ' - kinF^^) + h(cosF^^ - cosF^^)]i (3. 140) 


where 


n = ,/X (3.141) 

a 

The shadow delay time can then be calculated. We desire to know the eccentric 
longitude at which the thrusters are turned on, thus given the exit angle and the 
delay time, Kepler's equation must be solved for the eccentric longitude F^ at 
turn-on. If 


Co = ^ Fex - ksinFg^ -i hcosFe^ 


(3. 142) 


Then the required angle satisfies 

F = C + ksinF - hcosF^ (3.143) 

This transcendental equation can be solved by an iterative procedure. Let 

F = C (3.144) 

o o 

F . 1 = C + ksinF - h cos F (3.145) 

n+1 o n n 

2 2 

This iteration will converge since k +h <1. The iteration can be halted when 

(3. 146) 


I Vi 


- 1 2 

where c is some small number which we set to 10 . Thus the eccentric longi- 

tude at thruster turn-on is obtained. 


dF ex 

For the costate equation, Eq. (3. 137), — is needed at turn-on. 


dF 


and — — were already obtained (App. G). F_ satisfied Eq, (3. 143) so 
dz ^ 


dF, 


dC 


(3. 147) 


1 -kcosF^-hsinF^ d^ 
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(3. 148) 


dC 

o 

dz 


= ^ trp + n — — + (1-kcosF -hslnF ) —j 

T ^z ex ex dz 


Since 


trp = at +bT +c 




(2aT+b) 


From the definition of t in Eq. (3. 140) is obtained 


(3. 149) 


= 

dz 



kcosFg^ 


-hsinFg^ 


dF 


ex 

dz 


dF 

(l-kcosF.^-hsLnF^„) - — — (3.150) 
en en dz J n 


Inserting these partials into Eq. (3.147) yields 


dF 

dz 


dF dF 

- = — r(2aT+b+l)r^x — - (2aT +b)r ^ + (t„-(2aT+b)T ) 1 (3. 151 

r^. ^ dz dz ^ 


where 

r^ = 1 - kcosF^ - hsinF^ 

r = 1 - kcosF - hsinF 
en en en 

r^ = 1 “ kcosF - hsinF 
ex ex ex 


(3. 152) 

(3. 153) 

(3. 154) 


and since n is only a function of semimajor axis by Eq. (3. 141), the only nonzero 
element of is 


dn ^ _ 3 n 
^a "S' a 


(3. 155) 


3. 10 Power and Fluence: The Damage Function 

Power degrades as a function of accumulated fluence which is represented 
by a damage function, D(N). This damage function is given in Section 4. To 
calculate the costate equation, its partial derivative is needed. The damage 
equation is stated here, . 
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D(N) 


(3. 156) 


explc^ + C^{log^^Ny^\ 


where N is the accumulated equivalent 1 MEV electron fluence and and are 
constants for any particular spacecraft, but are functions of cell thickness and 
base resistivity as indicated in Section 4.3 in Eqs. (4.5) and (4.6). Then 




121og^Qe 

N 


11 . 


C2(logjoN)“D(N) 


(3. 157) 


This form of the damage function can cause numerical difficulties at the 
initial orbit if the initial N(t^) is set at too small a value and the integration time 
step is large since ^ looks like a gamma function with very high values for 
N < 10^^. Typically^ a value of for a low orbit would be around 10^^, so that 
it would be reasonable to pick a starting value greater than 10 . In the coding 

for the program, if the initial N is zero, then a revised initial N is calculated 
which is the "average" of the amount of fluence encountered on the first orbit, i.e. 


N(t ) . T 


N(t) 


(3. 158) 


where T is the orbital period. Facility also exists for inputing a nonzero initial 
value for N. 


3. 11 The Flux Equations 

An analytic equation for the flux as a function of spatial coordinates is 
constructed in Section 4. 2. It is repeated here, with slightly different notation. 


N = exp^ T A. 


(3. 159) 


where flux (N) is in equivalent 1 MEV particles, 

S = ln(R - 1) (3- 160) 

where R is the radial distance from the Ekrth's center in units of Earth radii, and 


5 2 

^i = rn*^i.j+5m 

j = l m=0 


LAT^’^SH™ 


(3. 161) 
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where LAT is geographic latitude in degrees, SH is shield thickness in mils and 
the K. are coefficients given in Tables 4-4, 4-5. The total flux is a sum of 

up to four equations like Eq. (3, 158), two for the front shielding, two for the 
back shielding, one of the two for electrons, one for protons. If the back shielding 
is infinite then the sum will have only two members. If the front and back shielding 
are equal only one equation like Eq. (3. 158) need be calculated for protons and one 
for electrons and then to get total flux, their sum would be doubled. In the 
remainder of this section we will consider only a single equation like Eq. (3. 158). 

Since the shield thickness does not vary for a particular spacecraft the sum 
Involving SH can be done once initially and Eq. (3. 161) can be rewritten 


A. 

L 


L K w-i’^ 1 = 1 5 

j=l 


(3. 162) 


where 

i=l 5 j = 1. . . . , 5 (3.163) 

Here w is latitude in radians, DTK is the degree to radian conversion factor. 

Thus Eq. (3. 159) can be calculated at any point on the spacecraft orbit as a function 
of latitude, radial distance, and the constant coefficients of Eq, (3. 163). The 
radial distance is just 

R = Jx\ + y\ (3.164) 

where Xj and were given in Eq. (3. 19) and (3. 20). The latitude is the arcsine 
of the angle between the vector pointing through the north pole and the spacecraft 
radius vector. This is given by 

w = sin~^f ^3 ^1 ^3 ^1 ^ (3.165) 

R 

where and g^ are the third components of the unit vectors given in Eq. (3.4) . 
They are functions only of p and q whereas X^, (and R) are functions of a, h 
and k. 
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In order to calculate the costate equations the partial of flux with respect 
to the orbital elements must be known. 




(3. 166) 


where 


az -^2 ‘J 


(3. 167) 


Let 


^ - i-2 

E. = E K. . wJ (j-1) 

' j = 2 'J 


(3. 168) 


Then rewriting Eq. (3. 166), 




aN 

a 


(3. 169) 


Now 


^ = _i_(x ^ + Y, 
(R-l)R^ ^ 


(3. 170) 


For h and k, and were given In Appendix B. For the semimajor axis, a, 

gz az 


^ 

aa (R-l)a 


(3. 171) 


Latitude Is not a function of a (the a's In the numerator and In the denominator 
cancel out In Eq. (3. 165))so that = 0. For the other elements 

a 1 r^^l ^^3 

d w - I f -L t <T _ + X - Y 

sz ^3 az 1 


Rcosw 




«i'3«i83> !llU 

T2 (.^1 aT" az - 

K 


(3. 172) 


47 



where and are functions only of hand k and the partials of and g were given 
in Appendix D. Here the necessary components are 



2 

T. 


1+P +q 


(qg3 


+ W 3 ) 


(3. 173) 


aq 


2p 

1+p^+q 


2 S; 


^ 2q ; 
I4p^+q^ 3 


agg 

aq 




(- pf -3 + 5.3) 


(3. 174) 


(3, 175) 


(3. 176) 


This completes the terms given in Eq. (3. 105). 


3. 12 Summary of State and Costate Equations 


In this section the full seven dimensional state and costate equations are 

summarized. We will assume that we have analytic functions for flux and the 

damage function without specifying their exact form. Only thrusting in an inverse 

square field will be considered. The effect of oblatene.ss is easily included by 

adding z j„ to the orbital element equations and - J 9 to the orbital element 

“ ^ dz 


adjoint equations, where these expressions were given in Section 3, 7 and App. F. 

We can consider a seven dimensional state composed of five orbital elements, 
mass, and accumulated particle flux. 


X = 


z 

m 

N 


(3. 177) 


Since m and N are varying slowly the first approximation of these quantities as 
well as the orbital elements can be considered. The unaveraged state equations 
are 


i "" mb yi-(z, F)u 


(3.178) 


m 


^(x, t, F.u) 
c 


(3. 179) 
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N 


f(z.F) 


(3. 180) 


and 


n D(N) 

P = S cos[0(z,u)] (3.181) 

R^a) 


For a particular orbit P is assumed to be zero for eccentric longitudes between 
Fg and Fj where is the eccentric longitude at entry into the earth's shadow 
and Fj^ is the eccentric longitude either at exit from the shadow or after a turn- 
on delay time after leaving the shadow. 

The unaveraged Hamiltonian is 


H = X'^x = 

where 


H 

= X z 

(3. 183) 

z 

— z — 

H 

= X^m 

(3. 184) 

m 

m 


■ V 

(3. 1 85) 


The averaged HamLltonian is 

H = *\ ^ H(x, F, F.'ulsfh, k. F)dF = ^ 

(3. 186) 

+ Xj^m(x.r, F,u)s(h, k, F)dF + ^ Xj^N(z, F)s(h, k. F)dF 


The state equations are 


j 

^Iz 


F 

i(5Z,F, F,u)s(h,k,F)dF = ^P(x,F. F.u)M(z. F)s(h. k. F)dF 

me 


(3. 187) 
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r^H _ r^2 




in(x, t, F, u)s(H, F, F)dF = - ^P(x, t, F,u)s(H, IJ, F)dF 


c ^F 


N - = \ N(z, F)s(H, Ic, F)dF 


The averaged Hamiltonian can be broken up into three parts 


H = H + H + 

z m N 


(3. 190) 


where 


H 

- vT- 

- A Z 

(3, 191) 

z 

— Z — 


H 

= X fn 

(3, 192) 

m 

m 


= V 

(3, 193) 


The costate derivatives are 


/xjM (X. t,F,u) . s(h,k.F) 
az -Jf^ az “ 


(Zzi^2£« t.F.u) + X^m(x,r. F,u) ^ (K, IE, F)1 dF 


cijr — T • — ^ ^ “i 

“ I ^2'-^ F.u)Js(E. 1E, Fg) 

F„‘ 


fiTT T— T . ^ ^ ^ “1 

~ t t, F., u) + X m(x, t, F, u) s(H, Ic, F- ) 

^z i-z m- -j 1 

■ 


5 "" ixN^(z.F)s{h,k,F)+Xj^N(z.F)^(H.lE.F)'ldF (3.194) 


•TT az 
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This can also be written 


F, 


2 ^ z A 

= - C I x^ -^(x.r.F.u)s(H.l^.F) + [H^(x.X2.*f^F.u) 






^s(H, Tt, F)dF 
dz 


(h (x.X ,t, F u) + t,F,u)s(H.Tc,F2) 


dz 


* 2 




dz . p 


Next 


^ 5N(z,F) 

\ 


r 


7T 




:(K,i;,F)+li<l.F)H(S:^}aF 

SZ 


X = - ^ = V \ xT P(x,r, F,u)M(z,F)s(h,F,F)dF 
1 


— m — — 2 ~z ' — ‘ 

Sm m c ' F 


but this just yields 




m 


Finally 




^ = - 2 r 


^2 aP(x,t, F.u) ,xJm(z.F)u 


aTT 


- - ^ls(h, k. F)dF 

r* 


m 


The dependence of P on N comes from the D(N) factor, the damage 
Thus simplifying 


aP(N) 

aN D(N) 


(3. 195) 
(3. 196) 

(3. 197) 

(3. 198) 
function 

( 3 , 199 ) 
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3. 13 Spacecraft Parameter Output 


It was desired to have the option of printing various spacecraft parameters 
at various points on an orbit. These parameters and the method of calculating 
them are discussed in this section. The parameters include: in-plane and out-of- 
plane thrust angles, yaw, pitch, and roll angles, the panel orientation angle, the 
sun incidence angles on the panel and the three contiguous sides of the center body. 

A nominal x-y-z body centered coordinate system was defined earlier and 
illustrated in Fig. 1-1. The yaw axis, e^, is pointed along the radius vector toward 
the earth, the pitch axis, e^ is normal to the orbit plane in a "downward" direction, 
the roll axis, e is in the orbit plane perpendicular to the radius vector in the 

^ A ^ 

forward hemisphere. In terms of the equinoctial orbital frame f , g, co. 


e 

— X 

If, 1a 

- f 4 g 

r — r 

(3. 200) 

e 

-y 

- w 

(3. 201) 

e = 
— z 

r — r E 

(3. 202) 


The thrust acceleration direction is available in the equinoctial coordinate 
system. Using the above transformation it can be obtained in the e^, e , e^ system. 
The in-plane angle, 6., and the out- of- plane angle, 6^, are related to the unit 
thrust acceleration^direction by 

u = cos cos - sin - cos 0^ sin e.e^ (3. 203) 


which Is illustrated in Fig. 3-7, where e 


- e and e. = - e 
-z -h -y 



Fig. 3-7 In-plane and Out-of-plane Angles 
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Then can be calculated from 


sin = - Uy . l0ol 


u 

^ - X 

cos 0. - 

' cos 9^ 


u 

• * - z 

sin 0. - 

cos 9 q 


(3. 204) 


(3. 205) 


(3. 206) 


If pitch and roll are constrained to be zero then 9 . = 0. 

The Euler angle order is yaw, pitch, and roll. The thrust acceleration 
direction as given in terms of yaw, >p^, and pitch, Q, by 

u = cos e cos + cos e sin i/iyCy - sin ee^ (3. 207) 

which is illustrated in Fig. 3-8. Note that the axes are not the same as in Fig. 3-7. 



Fig. 3-8 Yaw and Pitch Angles 
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since u Is known and q can be calculated from 


sm 9 = 




(3. 208) 


cos lb 


cos 9 


(3. 209) 


sin , 


cos 9 


(3. 210) 


For the pitch constrained case 9 = 0 and ip^ = 9^. 

For the constrained case roll is assumed to be zero. If there Is no roll 
constraint, then by a roll rotationand by a panel rotation about their longitudinal 
axis the panels can be faced directly at the sun to produce maximum power. If the 
roll is accomplished before the panel rotation, then the roll angle should be such 
that the panel axis is perpendicular to the sun line. The vector along the panel 
axis (also the pitch axis) after the yaw, pitch, and roll rotations should have an 
inner product with the vector pointing toward the sun of zero. Thus 


.„xyz.T T T r -,T 
[0] [a] 


1 

0 


= 0 


(3.211) 


where the yaw rotation matrix is 


[*y] 


COS Ip 

sin ip 
0 


sin ip 
cos Ip 
0 


(3.212) 


The pitch rotation matrix is 


[e] = 


and the roll rotation matrix is 


cos 0 

0 

- sin 0 

0 

1 

0 

sin e 

0 

cos 0 


(3.213) 


[^] = 


1 

0 

0 


0 

COS a 
- sin a 


0 

sin a 
cosa 


(3.214) 
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Since and 9 are already known let 




R' 

y 

R' 

z 


= [e]WylRr 


and then 


R ' cos a + R ' sin 0 = 0 

y Z 


(3. 215) 


(3. 216) 


Let 


tan 0 



(3. 217) 


with a 360° range. 

We assume the panel orientation angle is zero when the normal to the panel 
is pointed in the z direction. Since for the roll free case the panel axis is perpen- 
dicular to the sun direction after roll, the panel can be faced directly at the sun. 

If roll is constrained to be zero, the panel will be rotated to minimize the angle 
between the normal to the panels and the direction to the sun. Let 


R'- 


R" 

y 


[?il [el [0] 


(3.218) 


which is the sun direction in the rotated coordinate frame. Since the normal to the 
panel is nominally along the z direction it must be rotated by an angle given by 


cos r ' 


(3. 219) 


for the unconstrained case. The sun incidence angle will then be zero. For the 
constrained case, we want to maximize the cosine of the incidence angle. This is 
given by 


cos ip 


cos ^ R" + sin ^ R^ 


(3.220) 


which is maximized by setting 


cos ^ 


R 


ft 

z 



(3,221) 
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(3. 22 2) 


sin ^ 


R" 

X 




^ T 

R" 4 - R"^ 
X z 


and then 


^ (3. 223) 

The above three formulas are valid for the unconstrained case also since then 
R" = 0 and i = 0° . 

y P 

After the three Euler rotations the direction of the sun in the rotated system 
is just (R^» » ^ 2 ^' These are the cosines of the angle between the normal to 

the sides of the spacecraft and the spacecraft sun vector and thus yield the incidence 
angles on the spacecraft surfaces. 

Some additional quantities are also calculated such as the sun angle in x-y, 
i.e., the angle which the projection of the sun vector on the e^, e plane makes 
with the e^ axis and the thrust angle >li defined as tl> = - sun angle. This is the 

control angle used in Section 3.6. 
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SECTION 4 


RADIATION AND POWER LOSS MODEL 


4. 1 In troduction 

Solar cell degradation models are an integral part of low thrust solar electric 

spacecraft trajectory calculations. Eor these missions the thrust magnitude decays 

as the solar cell power is degraded by exposure of the cell to energetic trapped 

particles. Solar electric spacecraft trajectory computer programs include both 

codes that accurately simulate a mission and codes that optimize trajectories. An 

(5) 

example of the former Is the program SEOR which uses a detailed radiation 
degradation calculation based upon the models of Obenschain The code described 
in this report and Its predecessor in Ref, 2 are examples of trajectory optimiza- 
tion codes. These two types of codes place different requirements on the model used 
for radiation degradation calculations. For simulations the computational goal Is 
maximum realism so that an accurate assessment of the mission requirements 
can be made. For optimization codes the goal of maximum realism Is maintained, 
but the constraint of maintaining tractable run times imposes severe limitations on 
the degradation calculations. 

There are two basic types of degradation models which have been constructed: 
Those that recall basic NSSDC (National Space Science Data Center) subroutines for 
the magnetic field and spectral flux and determine the degradation integral at each 
trajectory sample point; and those which establish a data file of orbital average 
degradation. An example of the former Is the highly accurate and comprehensive 
degradation calculation used In SEOR^^^ The SEOR type procedure whereby the 
fluence is calculated during the trajectory integration via the calling of several 
NSSDC subroutines would be excessively time consuming when used in an optimiza- 
tion program. For an orbit of 18 sample points, and for 15 energy bands, the 
NSSDC subroutines INVARA and MODEL combined with a simple energy integration 
code require a computation time of 0.2 seconds /orbit on an IBM 360. A typical 
trajectory optimization run with SECKSPOT may comprise some 50 separate 
trajectories and iterate the solar cell degradation calculation 50, 000 times. There- 
fore several hours of machine time would be required If these calculations were 
made by the standard codes. Thus efficiency demands an approach for the degrada- 
tion model used in optimization codes different than that frequently used in simulations. 
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To achieve this efficiency most degradation models used in optimization 

programs are of the data file type. The model in the GEOTOP^^^ ^ program is of 

this type and is based on a table look up approach allowing only limited values of 

orbit parameters between which interpolation is used. Entries in the table are 

total IMEV equivalent fluence for one orbit, and the table arguments are apogee, 

perigee, and inclination. The maximum orbital inclination of the GEOTOP model 

is 30°. A simpler model than the GEOTOP model has been used in optimization 

(4) 

codes such as MOLTOP . This model is restricted to a single inclination and has 
no provision for changing cell characteristics. 

The degradation model presented in this paper extends the latitude range, 
provides a continuous and smooth representation of the flux field, and provides for 
changing the cell characteristics. And the deficiencies of the above data file models 
have been avoided without sacrificing computational speed; a typical trajectory 
fluence calculation requires only a fraction of a minute for thousands of sample 
points. Further, the model is extremely simple to code, with all required informa- 
tion presented in this section. It is based upon two analytic expressions: D(N, BR, TH) 
which describes solar cell power degradation as a function IMEV equivalent fluence 
(N), cell base resistivity (BR) and thickness (TH), and f(R, L, SH) which describes 
a spatial field of IMEV equivalent electron flux. The parameter R is distance from 
the earth's centroid, L is latitude, and SH is cell shield thickness. When the flux 
f is time integrated over a sequence of trajectory time steps the result is a fluence 
value N which is entered into the degradation equation D(N, BR, TH) to obtain the 
fractional power loss. By constructing the model in an analytic format it is possible 
to make computation time minimal, and also to allow analytic differentiation of the 
flux and power loss for computation of adjoint variables for optimization studies in 
codes like SECKSPOT. The analytical expressions for D and f contain a number 
of coefficients derived externally from the trajectory optimization code using the 
NSSDC codes. A major advantage of this approach is that the coefficients can be 
altered to reflect new cell damage data or radiation belt data with minimal effort. 
IMEV flux values are averaged over one day in the model. Although this is not a 
requirement for the construction of this type of model it eliminates the requirement 
for integration time steps that are a small fraction of a day. The flux field is 
defined in geocentric coordinates to eliminate the need for a series of transformations 
from geomagnetic to geocentric coordinates. A model dividing the problem into the 
same two sub-models was previously constructed by Ives^^^^, however, the Ives 
model was tabular in format, limited to equatorial orbits, and did not contain cell 
characteristics as parameters. 
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4. 2 IMEV Equivalent Electron Flux Field Model 

Construction of the IMEV flux analytical model is a two step procedure. 

In the first step a space of discrete values of IMEV flux is generated from particle 
spectral flux data, and from data concerning the conversion of spectral flux to 
IMEV flux. The second step consists of fitting the array of discrete points with 
analytic expressions by means of a multivariate regression analysis. 

Generation of the IMEV equivalent electron field points is accomplished by 
a series of integrations. The complete procedure is illustrated in Fig. 4-1. A 
symbolic guide to the flow diagram is given below in Table 4- l.The process begins 
by calling the NSSDC code INVARA^^^^ which converts a set of geocentric coordinates 
to the magnetic coordinates B and L. B is the magnetic field intensity, and L is the 
magnetic invariant associated with particle momemtum along the magnetic field 
lines. Seven choices of a model for the magnetic field are available within the code 
INVARA. For the work here the IGRF (International Geophysical Reference Field) 
option has been selected. This option was recommended by E. G. Stassinopoulas 
(in a private communication) of NSSOC as an often used reference model. The time 
selected for the field was 1978. 
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Loop on Coordinates 























Table 4-1 Legend for Figure 4-1 



Indicates code, activity 


Input Function 


Branch point 


Definitions : K 

6 

B 

L 

EE, I P 

F 

SH 

' 1 ME\ 


Results 

Radius, i.e., dist. from center of earth 

Latitude 

Longitude 

Magnetic field intensity 
Magnetic invariant 
Electron and proton energy array • 
Particle flux (spectral) 

Shield thickness 

Equivalent 1 MEV particle flux 


Codes: INVAI'A 

MODEL 
MAIN 
DAM E 
DAM 

REGR'^SS 


From NSSDC, magnetic field code 

From NSSDC, Hux code 

Calling code, input output, etc. 

Electron damage subroutine 

Proton damage subroutine 

IBM code for multivariate regression 


61 





The magnetic coordinates B and L are entered into the subroutine MODEL 
which computes proton and electron flux within a selected energy interval. This is 
done for the entire range of important energies. MODEL^^^^ is a standard particle 
flux code produced by NSSDC, and the fluxes it produces are based upon smoothed 
satellite data. The data in MODEL is in block form, and interpolation is used 
between defined points. MODEL data is based upon satellite measurements made 
in the period 1964-1967, It represents a solar cycle maximum. The energy 
intervals selected for flux evaluation are dictated in part by the MODEL code. The 
code allows energy increments which vary over the total energy range. Practical 
upper limits on energy are set by flux values which are less than unity. Lower 
energy limits are set both by MODEL and by shielding considerations. For protons, 
the lower energy cut-off is related to shield thickness (SH) by the relationship 

EP (MEV) = 1.53 

f35) 

which was obtained by a fit to the data of Rasmussen . For both electrons and 
protons fifteen energy steps have been used. The proton energy range is 0 to 40 
MEV with the lower end adjusted by the above formula. For electrons the range 
used was 0. 2 to 4, 8 MEV. 

Spectral flux values produced by subroutine MODEL are averaged over 
longitude to give a one day mean flux. The aim of this averaging is to avoid keeping 
time of day in the trajectory optimization code. Retaining time would require 
additional coordinate transformations and would add another dimension of complexity 
to the field fitting problem. Variation of the flux field with longitude is considerable. 
An example of the range of flux values is shown below in Table 4-2 for protons in 
the energy range 4-5 MEV. 

Table 4-2 Longitudinal Variation of Proton Flux 
(par tide s/cm2 /sec) 


Longitude 

0° 

60° 

120° 

Flux 

. 135E06 

. 540E06 

. 871E06 

Longitude 

180° 

240° 

300° 

Flux 

.393E06 

. 794E04 

.449E03 


These values have been obtained near the proton peak at 2. 2 earth radii, and at a 
North latitude of 21. 5® The average value is . 330E06. Variations with longitude 
are caused by two effects, namely the non-dipole asymmetry in the magnetic field, 
and the non-alignment of the magnetic pole with the geodetic pole. Generally 
asymmetry effects dominate. 
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These values have been obtained near the proton peak at 2. 2 earth radii, and at a 
North latitude of 21. 5°. The average value is . 330E06, Variations with longitude 
are caused by two effects, namely the non-dipole asymmetry in the magnetic field, 
and the non-alignment of the magnetic pole with the geodetic pole. Generally 
asymmetry effects dominate. 

Orbital fluence calculations are little affected by the longitudinal averaging 
since the rotation of the earth and the consequent rotation of the radiation belts 
does effectively average the fluence. 

The next step In the field model construction is an Integration which converts 
the azimuth-averaged spectral flux data to an equivalent IMEV electron flux. This 
is accomplished with a weighting function in a Fredholm Integral. Using the 
Figure 4-1 notation, the integral for electrons) is 

F'(SH) = ^ F(EE) DAME(EE, SH) dEE 

where EE is the electron energy and SH is the shield thickness, DAME is the 
electron weight function and F' is the spectral, longitude averaged flux. The 
function DAME was constructed from data taken from the curves of Rasmussen 
however because of the inaccuracy of this type of data transfer the data had to be 
restructured. It can be argued by analogy with optical extinction processes that 
the relationship between shield thickness and the weight function should be 
exponential, i. e., 

log.nDAME = - m(EE)SH + log.^DAMEl 

<=10 lu 'SH=0 

This function was found to fit the data very closely and thus was used for the 
restructuring. The functional form for m(EE) was found to be 

log^Qm(EE) = m' . EE + C 


where m ' and C were obtained by regression. In terms of energy and shield 
thickness, the function DAME can thus be written 


DAME 


B(EE) 10 


-dO”^ 


. EE+Cjgjj 


The coefficient B(EE) is the appropriate DAME function for SH=0, Again, by a 
regression analysis, B(EE) was found to be (to ± 3%) 
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The equivalent IMEV flux for protons is computed by the Fredholm integral 
F ' (SH) = \ DAMP(EP. SH) F(EP)dEP 

where EP is proton energy. The weighting function for protons (DAMP) was also 

(35) 

obtained from the curves of Rasmussen . In addition to the cut-off energy 
described above, the function DAMP is defined by 


DAMP = 


jQ-O. 721og^QEP + 4. 13 

EPq s EP < 11 MEV 

2200 

11 < EP < 46 MEV 

llog2QEP“^5. 26 

46 < EP MEV 


Computation of the weight integrals completes the process of generating raw 
data for the field model. The raw data is a collection of IMEV equivalent electron 
values defined at discrete spatial points with solar ceil shield thickness as a parameter. 
Graphical examples of the flux field data are shown in Figures 4-2 and 4-3. Points are 
generated at 10 ^ latitude increments for the radii listed in Table 4-3, IMEV flux 
was computed for four values of shield thickness, namely 3, 10, 20, 30 mils. The 
various combinations of latitude, radii, and shield thickness yield approximately 
2500 flux values. 



Fig. 4-2 Electron Contribution to Total Flux vs. 
Radius, SH=10 mils 
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The second part of the flux generating process consists of fitting these discrete 
points with analytic functions. A multivariate regression (subroutine STEPR)^^^^ was 
used as the vehicle for the data fitting, and a two step process was used for protons 
and electrons separately. In the first step the radial flux points were fit for each 
latitude and shield thickness. Since the electron data has two peaks at the inner 
and outer Van Allen belts, quartic functions of radius were used. To reduce the 
dynamic range of the flux, In(flux) was used as the dependent variable. Finally 
to force the flux to zero as radius approaches unity, the radial variable S = In(R-l) 
was used. Thus the radial functions are of the form 

4 

ln(F) = L A.S‘ (4.1) 

i=0 '■ 

where the coefficients A^ are functions of latitude and shield thickness. These 
coefficients were then fit with a function of latitude and shield thickness, vis., 

4 2 

A: = L L (B..ShJ)C.. LAT^ (4.2) 

^ k=0 j = 0 

Substituting (4.2) in (4. 1) for A. and taking the exponent thus yields a single analytic 
expression for the IMEV equivalent electron flux as a function of latitude, radial 
distance and shield thickness. There are four separate Eq. 4. I's, one for protons 
and one for electrons for both the front and the back shielding, and the fluxes from 
these are added to obtain the complete flux. The coefficients used are tabulated in 
Tables 4-4 and 4-5, however the regression analysis yields directly products of 
B and C. and therefore the products are listed instead of B and C. We redefine A. 

i 

as follows: 

5 2 

A: = L E K ... LAT->"-'sh”" (4.3) 

^ j = l m=0 

This field model is not completely general. In the latitude region 50® to 
60® it is valid only for R > 1. 5 earth radii. In the region 60® to 75® A is 
restricted to R > 2. 5. The restrictions stem from function fitting problems and 
not from the data. Hopefully they can be removed with improved function choices. 
Above a latitude of 75® the field is set to zero in accord with the zero flux output 
from subroutine MODEL. 
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Table 4-4 Expansion Coefficients for IMEV Electron Contribution to 
Total Flux (Second line starts with K(6) etc. ) 

i=0 0. 1 0 ? 0. -0. 77‘^)t>;'BH7F-03 -n . ] 3 1 04^3F-0^ o 

-0. 07F-01 0. ^^OO-^^nSF-OA 0 . 1 *>8630 1 AF-0^. -O. 1 7/t?3() ) F>F-07 -0 . 3 8 F -08 

0. 131 1 7390F-03 -0. 1 7f'6?F84F-0b 0. 1 ?3F>5317F-00 0. 3R83F078F-00 -0 . 3 ? F ^7 8 ?7 F - 1 0 

i= 1 -0 ni -0.4b875O44E-0? 0 . 3b49 53?0E-0? 0 . 3350P00 1 F-00 -0 .0? 2 1 3 1 30 F-08 

0.927980b8E-n? -0 . 0 7 3 0 ] 1 89 F -03 0. 1 8 1 91 829F-04 0 . 1 89 1 8 797 E -08 -0 . 88 3 8 0078 F -0 w 

-0.29772083F-03 0 . ? 1 207 1 78 E-09 -0.49198288F-O7 -0.80891998F-OH 0 . 38 8 ?8 3 98 F - 1 0 

i=2 0.22939828F 01 -0 . 807 9 29 8 9 F -0 3 -O . | 3 2389 9 8F -02 -O . 8 1 278 1 1 3F -08 -0 . 289 87 2 7 ? F-0 7 

0.28938293F-01 -0 . 8 79 3 1 7 9 2 F -09 -0 . t 8888 1 98 F-f)9 0 . 9 8 8 7 3998F -08 0.188270888-08 

0. 27908 ] 98F-03 0.937809378-08 -0 . 3 1 999 8 1 1 F -08 -O. H7088 1 89F -09 0 . 8 8 7 3<-t8 1 8 F - 1 0 

i= 3 0.88829821F 00 0. 12882998F-01 -0 . 11 7 8 7889 F-02 -0 . 908900H9F-0 8 0 . 1 9 28 2 1 79 F - 08 

-0. 80808879F-02 0 . 9 0 1 8 39 2 8 E -03 -0 . 87 283309F -08 -0 . ] 1 1 9888 3F -08 0 . 1 8 3 ] 8 08 9 F -08 

0. 1 3798090E-03 -0. 1 2887829E-09 -0. 1 1 8793 1 7F-07 0.38839289E-08 -0. 1 2 183820F-1 0 

i= 4 -0.97988923F 00 - 0 . 8 3 ] 2 22 8 7 F - 02 0 . 9 3088«8 1 F -03 O. 1 87 7 1 09 1 F -o8 -'). 70028888 F-07 

-0.89279007F-f)2 0 . ] 2098987F -09 O .] 998988 7 F-08 0 . 9 1 1 99928 F -08 -0 . 99 8 399 7 3 F - 1 (' 

-o. 122 1 8909F-03 -0. 1 338 1 789F-08 O. 1 9723730F-08 O. 1 70998 1 OF -09 -0. 3 1 3 802 28 F - 1 0 


Table 4-5 Expansion Coefficients for Proton Contribution to 
Total Flux 


i=0 


i= 1 


0.2A5505fi3E 02 
-0.39871807E 00 
0.89885H09E-02 


-0.R3898939E-02 
0.201 H99H3F-02 

-0.92885 19 7F^-09 


-0.19293925E-02 
-0.7803177BE-09 
0. 198 28506E-0 5 


0.2801R12AE-05 
-0. 14334919E-05 
0.30082988E-07 


-0.37385fl82E-08 

-0.93087397F-07 

0.71108997F-09 


0.13719120E 01 
-0.3R838380E 00 
0.78533999E-02 


-0.29988785E-02 

0.20878092F-02 

-0.20R799H9E-09 


-0.4892R572E-03 
-0. 10R578R5E-03 
0.539fl93A8E-05 


0.112A7R25E-09 
-0. 19893597E-05 
-0.12A38271E-07 


-0.1009388RF-08 
-0.79271 8R2F-0B 
-0.27379979F-0R 


i=2 -O.38R0R71OE 01 0. 57829099F-0 1 0.87337159E-09 -0.9899 1385E-09 -0.R393R189E-08 

-0.1R90791)E 00 -0.7277622RE-02 -0. 38R22810E-0A 0. 78R09995E-05 0.9 2328708E-07 

0. 389R9293E-02 0. 29P8H8 1 3F-03 0 . 8885778 9E- 05 -0.31987960E-08 -0. R590838RF-0R 

i= 3 -0.21933759E 01 0. R02309 1 5E-01 0.3929R80RE-03 -0.800593996-09 -0. 10387 207F-08 

0.1R117610E-01 -0.70728998F-02 0.358808616-04 0.48947443E-05 0.22437872F-08 

-0.23797412E-03 0. 12794089E-03 -0.2B894892E-05 -0. 5 1 577029E-07 0. 1 1804R97E-0H 

1=4 -0.44003046E 00 -0. 144H5885E-01 -0, 1 2907069E-02 -0.22752400E-04 0. 10778589E-06 

0.85140910E-02 -0. H78888486-05 0. 988541706-04 0. 2679 105BE-05 0 . 19 1 53641 F-08 

-0.20893010E-03 0.225H6333E-04 -0. 148 13049E-05 -0.972171966-07 -0.115993226-08 




67 



Over the defined region of R, LAT, SH space the field accuracy is as follows: 
For protons the mean deviation of the value predicted by Eq. (4-3) from the MODEL 
value at the defined R,LAT points is 5. 6 %. The average of the absolute value of 
the deviations is 28%, For electrons the mean deviation Ls 1. 2% and the average 
absolute deviation is 25%, 

4, 3 Power Loss Model 

The second part of the two-part analytical model is an expression relating 
IMEV fluence to fractional power loss. The data for this expression is provided 
by the curves of Carter and Tada , and it can be fit closely by an expression of 
the form 


InD = + C 2 [logjQSN]^^ (4.4) 

Here D is the fractional power loss after the cells have been subjected to fluence 
N ( 5 ^ N in units of IMEV equivalent electrons/cm /sec and is summed over front 
and back shielding and electrons and protons). Constants C^ and C^ are functions 
of cell thickness, and a simple quadratic has been found to adequately represent 
this variable. Although variations in cell base resistivity have a relatively small 
effect on the degradation function, the effect is represented in the Carter and Tada^^^^ 
data. Two sets of degradation curves are presented, one for base resistivity 1 to 
3,Qcm, the other for 7 to ISQcm. We therefore provide two sets of functions for 
Cj and C 2 » TH is the cell thickness in mils. 


l-3Qcm 


7-130cm 


C 


2 


[0. 22647 + 0. 05217 TH - 0. 00443 (TH)^]0. 1 
[+ 05151+. 03641 TH-. 00144 (TH)^](-10'^'^) 
[0.04914+0.07056 TH - 0. 00435 (TH^] 0.1 
[-0. 05118 + 0. 04517 TH-0. 00156 (TH)^](-10“^^) 


(4.5) 


(4.6) 


The C values are defined for cell thickness in the range 4 to 14 mil, and D is valid 
for fluence values above 10 ^^ (> 0 . 1 % degradation). 
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4,4 Conclusions 


Since trajectory optimization normally involves the evaluation of several 
quantities which are functions of spatial coordinates at a sequence of time steps, 
it is a simple matter to apply the analytic degradation model described here. At 
each sample point in the trajectory the flux function (4. 3) is evaluated. Multiplication 
of the resultant flux by the sample time step yields an incremental fluence, and 
these are simply summed up to a given sample time to provide a total fluence which 
is used in (4.4) to find the power loss. Since this process involves the evaluation 
of only two analytic expressions it requires very little computation time. Modeling 
of the IMEV flux field as a separate entity allows simple consideration of both 
front and back shielding. The various coefficients in the analytic expressions 
relate to specific cell damage data, however having established the general 
analytical characteristics of the model, it is a simple matter to update the co- 
efficients using the latest cell damage data. 
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SECTION 5 


NUMERICAL RESULTS 


Two sets of cases will be discussed Ln this section. The first, which assumes 

an inverse square gravity field and thrusting with no oblateness, degradation, or 

(25) 

shadowing was also reported previously . This set includes both attitude 
constrained and unconstrained examples. The second set includes the new degrada- 
tion model reported in Section 4 as well as other perturbations. All examples in 
this set include attitude constraints. All cases were run on an IBM 360/75 computer. 
The code is in Fortran IV using double precision. 

5, 1 Examples Without Power Degradation 

Several SERT-C type cases were run. All of these assume only the inverse 
square gravitational field with or without attitude constraints. No oblateness, 
shadowing or power degradation due to radiation are assumed. Comparisons are 
made between the unconstrained case and attitude constrained cases while varying 
launch date and time. A particular constrained example is looked at in greater 
detail. 

Table 5-1 summarizes the characteristics of the SERT-C type cases. The 
final orbit is geosynchronous. These cases were run with a 10 day time step and 
averaged around an orbit using either two 4 point or two 8 point quadratures. The 
variation in resulting flight times and ^V's from using sets of 4 or 8 point 
quadratures was less than , 5%. 


Table 5-1 SERT-C Example Data 


Initial a 

= 

9528. 16 km 

Initial e 

= 

0. 

Initial i 

= 

28.3® 

Initial mass 

= 

849. 6 kg 

Maximum power 

= 

4. 828 kw 

Specific Impulse 

= 

2900 sec 

Final a 


42164. km 

Final e 

= 

0 

Final i 


0® 
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The effect of varying the time of launch during a day for the constrained case 
can be shown by varying the Initial longitude of the ascending node. Changing the 
launch time affects the sun-spacecraft geometry and therefore the resulting 
trajectory. This effect Is Illustrated in Fig. 5-1 for a launch date of March 21. 5, 1980 
The flight time (t^) and are plotted versus initial longitude of ascending node. 

The dotted line Indicates the unconstrained value. The flight time varies by about 
18 days with a minimum of 130 days around p = 135° compared to the unconstrained 
value of 124 days. The range of increased flight time over the unconstrained case 
is from about 5 to 20%. The unconstrained is 4.65 km/sec and the constrained 
cases vary from less than 1% to about 12% greater. In the minimum case about 
15% of the initial mass was consumed. This was only about . 6 kg more than the 
unconstrained case. 



Fig. 5-1 Flight Time and Versus Nodal Angle 



Four cases of the attitude constrained example were run varying the date of 
launch. The initial longitude of ascending node was 0°, The percent increase in 
flight time and over the unconstrained case is shown in Fig. 5-2. For any date 
of launch a curve such as shown in Fig. 5-1 would result. Although an absolute 
minimum time would occur on a particular day, it is expected that by varying the 
time of launch during the day that a nearly minimum flight time would be found. 

In summary, for this example, for appropriate selection of launch date and time the 
attitude constrained case results in as little as a 1% increase in and a 5% increase 
in flight time over the nonconstrained case. 



21.5 21.5 21.5 21.5 

Fig. 5-2 Flight Time and for Various Launch Dates 

The o = 135° case which results in the lowest flight time and will be 
looked at in greater detail, including further comparisons with the unconstrained 
case. The semimajor axis and inclination histories for these two cases are 
plotted in Fig. 5-3. The maximum yaw angle is plotted in Fig. 5-4 for the two 
cases. The yaw deviates from zero by plus and minus approximately this same 
amount over one orbit. For the unconstrained case the maximum yaw occurs when 
the spacecraft is at the line of nodes. As the spacecraft nears geosynchronous 
orbit the maximum yaw increases to over 90° for the constrained case, but to 
somewhat less, 73°, for the unconstrained case. Since only discrete points every 
15° were printed out, the actual maximum may have been slightly higher than 
plotted. For this circle to circle transfer the pitch angle for the unconstrained 
case is always zero. 
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Fig. 5-3 Semimajor Axis and Inclination Histories 



Fig. 5-4 Maximum Yaw Histories 
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Fig. 5-5 shows yaw plotted as a function of eccentric longitude for the initial 
and the final orbits. Because of the coupling with power, the maximum does not 
occur at exactly the line of nodes (which would be 135° and 315°) as in the uncon- 
strained case. There are two jumps in yaw of about 60° for the final orbit. This 
will be discussed further. 



Fig. 5-5 Yaw for Initial and Final Orbits 


Fig. 5-6 shows a plot of the ratio of power to maximum power versus 
eccentric longitude for the initial and final orbits. For the initial orbit the power 
is near maximum for nearly the entire orbit, dropping to a minimum of . 91 but 
with an average of . 98. For the final orbit the minimum reaches . 707 which is the 
absolute minimum possible and which occurs at a jump point of the control. Even 
so, the average power ratio is over . 9. The two other relative minimums for each 
orbit occur near to the minimum and maximum of the primer vector angle (a) (and 
also the thrust angle, ip). 



Fig. 5-6 Ratio of Power to Maximum Power for 
Initial and Final Orbits 
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In Fig. 5-7 are plotted the primer vector angle (a), the thrust angle (}b), 
and the sun angle (^) for the initial orbit. These angles are defined in the “ £2 
- e^ system defined in Section 3.5 and were illustrated in Fig, 3-2 and Fig, 3-4. 
The sun angle is the angle between the radius vector and the line to the sun, and 
can vary by at most 180° depending on the orientation of the orbit. For the initial 
orbit the variation is approximately 140° , The thrust angle is a function of ci and 
p and is equal to ^ at “0° and 180°. For this orbit = 90° when = 90°. Since 
g 157° and 23° for the two positions at which a - 90° there is no jump in ^ (for a 
jump 45° < p < 135°). 



Fig. 5-7 Primer Angle (a). Thrust Angle (j/)) 
and Sun Angle ip) for Initial Orbit 


The primer vector angle (^), thrust angle (^), and sun angle (g) are plotted 
for the final orbit in Fig. 5-8. Obviously for this orbit there is a jump in when 
a - 90°. Note that at this point g 124° for the first jump and g 56° for the 
second. For these values of g the jump in j/) is approximately 60°. It is interesting 
to note that after the first jump decreases slightly even as cy increases (0 in- 
creasing as a decreases after the second jump) as the changing p allows a closer 
alignment with the primer vector at less power loss penalty. Jumps in j/> occur 
during about the last 20 days of the trajectory. 
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Fig. 5-8 Primer Angle (a). Thrust Angle (ip) 
and Sun Angle (p) for Final Orbit 

The panel orientation angle (defined in Eq. (3. 85) is plotted in Fig. 5-9 
for the initial and final orbits. In both cases there is a rotation through a full 360°. 
For the final orbit there are discontinuities at the point where there are discontinuous 
changes in the thrust direction. Both jumps are about 80° (note panel angles of 180° 
and -180° are equivalent). 
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Fig. 5-9 Panel Orientation Angle for 
Initial and Final Orbits 


In summary, for the attitude constrained and unconstrained cases compared 
here, flight time and trajectory histories differed only somewhat. The main 
differences are the discontinuities in thrust direction and in solar panel orientation 
which occur during the last 20 days of the constrained case. Average power levels 
were at most 10% less than the unconstrained case. Since these were circle to 
circle transfers the effect of a zero pitch constraint Is not illustrated. 

5. 2 Examples With Power Degradation 

Several cases were run to illustrate the use of the final version of the code. 
Most of the cases were of a SERT-C type, although some GEOSEPS-type missions 
and some cases which originated in an elliptical orbit were run. All of these cases 
included oblateness and the new power degradation model. Unless otherwise noted, 
all assumed that pitch and roll were constrained to be zero. The constants in the 
coefficients of the power loss function for these runs were slightly different than 
those given in Eqs. 4. 5 and 4.6. 

The SERT C cases are discussed first. Table 5-2 lists conditions on the 
initial and final orbits and the spacecraft. 
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Table 5-2 SERT-C Example Data 


Initial a 

= 

9528 km 

Initial e 

= 

0 

Initial i 

= 

28. 3® 

Mass 

= 

850 kg 

Maximum initial jet power 

= 

4. 828 kw 

Specific Impulse 

= 

2900 sec 

Final a 

= 

42104 km 

Final e 

= 

0 

Final i 

= 

0° 


Cases with other final orbits were also run. 

i<‘or tne cases discussed in this section a four-point quadrature was used 
between the points where the primer vector angle, a , crossed + 90° (see 
Section 3.5). For the initially circular cases there were two intervals per orbit. 

A ten day time step was generally used, although some cases had a 15 day time 
step. A Newton-Raphson iterator and a four step Runge-Kutta integrator without 
accuracy controls were used. 

The characteristics of one set of runs are listed in Table 5-3. Launch 
date, launch time (determined by the initial longitude of ascending node, p ) and 
solar array parameters were varied. Shadowing was not included. A 10 ohm-cm 
solar cell base resistivity was assumed. If additional values for launch date and 
time were run then a curve such as in Fig. 5-1 could be drawn. For the limited 
number of runs shown in Table 5-3, the variations in transfer time and are 
small. The power degradation which is near 56% increases transfer time 
significantly compared to the constant power case discussed previously. Without 
infinite back shielding transfer time is increased by 14% as the degradation factor 
decreases an additional . 08. Increasing the cell thickness changes the shape of 
the degradation versus flux curve. The increase from 6 to 10 mils caused a final 
damage factor decrease from .56 to .41 and a 30% increase in transfer time. 

The case shown in the first column will be illustrated in more detail. 

Fig. 5-10 is a time history plot of semimajor axis, inclination, and the longitude 
of the ascending node. The change in the latter is largely due to oblateness which 
is strongest at lower altitudes and high inclinations. The semimajor axis actually 
overshoots the final value somewhat. This overshoot comes about because initially 
greater effort is put into reaching higher altitudes at higher inclinations where the 
flux is less. As the higher altitudes are reached, effort is expended to reduce 
inclination. 
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Fig. 5-10 Semimajor Axis, Inciination, B,nd 
Longitude of Node Histories 


Table 5-3 SERT C Example Results 
(10 ohm-cm base resistivity) 


Launch date 

Mar 21. 5 

Mar 21. 5 

Dec 21.5 

Dec 21.5 

Mar 21.5 

Mar 21. 

Longitude of node 

Qo 

45^ 

0^ 

45^ 

0^ 

0^ 

Cell Thickness (mils) 

6 ’ 

S 

6 

6 

6 

10 

Front Shield Thickness (mils) * 

6 

6 

6 

6 

6 

6 

Back Shield Thickness (mils) 

infinite 

infinite 

infinite 

infinite 

6 

infinite 

(km/sec) 

4. 89 

4. 92 

4. 82 

4. 91 

4. 85 

4. 83 

Transfer Time (days) 

236. 

236. 

231. 

237. 

269. 

308. 

nif/mo 

. 842 

.841 

. 844 

. 841 

. 843 

. 844 

Damage factor 

. 563 

. 571 

. 564 

. 569 

. 482 

. 412 

Initial Costate 

3657. 3 

8366.5 

4331. 1 

4797. 6 

2545. 9 

6074. 5 


121. 6 

286.5 

97. 6 

350. 0 

170. 7 

238. 6 


298.0 

86. 6 

263. 7 

71. 7 

441. 1 

615. 6 


-6072.8 

-25714. 

-6300. 6 

-26335. 

-8327. 

-8408. 9 


-35259. 

-24085. 

-33507. 

-19704. 

-31798. 

-29438. 


-33428. 

-37246. 

-31866. 

-31619. 

-35878. 

-46571 


-106. 87 

-121. 89 

-87.53 

-89. 73 

-52.31 

-174, 51 
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Fig. 5-11 shows the yaw angle for the initial and final orbits. For the initial 
orbit the variation is fairly smooth, with a range of about 60° either side of the 
orbit plane. The final orbit contains two jumps in the yaw angle; one at 1.3® and 
one at 180.3®. On each orbit of this trajectory there were two points at which the 
primer vector angle was ± 90® (see Section 3.5), but the sun angle was between 45® 
and 135® (thus causing a jump in yaw) for the first time at about 105 days. The 
final yaw angle has a range of 125° either side of the orbit plane. The thrusting at 
yaw magnitudes greater than 90® causes a decrease in semimajor axis near the end 
of the trajectory as illustrated in Fig. 5-10. 



Fig. 5-11 Yaw for the Initial and Final Orbits 


Fig. 5-12 shows the power variations caused by the inability of the panels 
to directly face the sun at all times on the initial and final orbits, where the 
maximum is unity. For the initial orbit power decreases to . 90 of its maximum, 
and the average is . 954, For the final orbit power decreases to . 707 of the maxi- 
mum at the jump points and the average is 0, 88, 



'#Lg. 5-12 Power Variation for the Initial and Final Orbits 
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In Fig. 5-13 are plotted the primer vector angle, the sun angle, g, and 
the thrust angle ij) for the final orbit (see Section 3, 5, especially Fig. 3-3 and 3-4). 
The primer vector angle and thrust angle coincide at 0^ or 180'^. When ot crosses 
-90° there are jumps of about 36° when g = 48° and g = 132° . 



Fig. 5-13 Primer Vector Angle (a), Sun Angle (fl), 
and Thrust Angle (^j) for the Final Orbit 


The equivalent of 1 MEV fluence is plotted in Fig. 5-14 for three cases. 
Power versus time is shown for the same cases in Fig. 5-15. The case with infinite 
back shielding and a cell thickness of 6 mils (first column of Table 5-3) has the 
lowest accumulated fluence and the least degradation (to .56) and the lowest transfer 
time. With a six mil back shield the total fluence is doubled and power degrades to 
.48. With infinite back shielding and a cell thickness of 10 mils (6th column of 
Table 5-3) the fluence curve coincides with the 6 mil case through 60 days, then 
increases somewhat more since at that time in the flight it is at lower altitudes and 
thus in a stronger radiation field. Its power curve in Fig. 5-15 falls the fastest to 
a final value of ,41. Note that most of the fluence is encountered during the first 
80 days or so and that power drops off sharply initially and reaches an essentially 
constant value after 80 or 100 days. 
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Jet Power (kW) Equivalent Fluence Particles) 
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6 mil front and back shield, 
6 mil cell thickness 



Time (Days) 


Fig. 5-14 Fluence versus Time 



Fig, 5-15 Power Degradation Histories 


82 





In order to further illustrate the effect of various array parameters, several 
cases of considerably shorter transfer time are listed in Table 5-4. These have the 
same initial data as shown in Table 5-3, but terminate at a = 12800 km, i - 25° . 

The third, fourth, and fifth columns correspond to shorter versions of the cases 
shown in the first, fifth, and sixth columns of Table 5-3. 


Table 5-4 Short SERT-C Example Results 


Launch Date 

Mar 21. 5 

Mar 21.5 

Mar 21.5 

Mar 21.5 

Longitude of node 

0 ^ 

0° 

0<^ 

0° 

Cell Thickness (mils) 

6 

6 

6 

6 

Front Shield Thickness (mils) 

20 

6 

6 

6 

Back Shield Thickness (mils) 

infinite 

infinite 

infinite 

6 

Base Resistivity (ohm-cm) 

10 

2 

10 

10 

iV 

1. 03 

1. 03 

1.03 

1. 03 

Transfer Time 

35.4 

46. 5 

41. 9 

46. 1 

Mf/M^ 

. 965 

.964 

. 964 

. 904 

Damage factor 

. 797 

.516 

. 611 

. 533 


Mar 21. 5 
0 ° 

10 

6 

6 

10 

1.03 
50 . 6 
. 964 
.462 


The case shown in the first column has a very thick front array shield of 20 mils. 

The degradation is about one half of that for the 6 mil shield listed in the third 
column. The ^V's are the same but flight time decreases by 15%. The case in the 
second column assumes a base resistivity of 2 ohm-om rather than 10 ohm-cm which 
was used for all other cases. This assumption affects the power versus fluence 
curve resulting in a greater degradation (an extra 10%) when compared to column 3. 


To further illustrate these five cases fluence is plotted in Fig. 5-16 and 
power in Fig. 5-17. The accumulated fluence does not vary with changes in cell 
thickness or base resistivity on the scale shown by the plot. After 40 days, the 
case with infinite back shielding and 20 mil front shielding has about 10% as much 
fluence as the case with 6 mil front shielding and 5% as much fluence as the 6 mil 
front and 6 mil back shielded case. Power variations diverge for the five cases, 
some of which reach near-steady state values in 40 days. However, the 20 mil 
shielded case is at a much higher power level. The 10 mil cell thickness case 
results in the greatest degradation. 
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Jet Power (kW) Equivalent Fluence (10 ° Particles) 


A, 20 mil front, infinite back 
shield, 6 mil cell thickness, 
10 ohm-cm base resistivity 


B. 6 mil front, infinite back 
shield, 6 mil cell thickness, 
2 ohm-cm base resistivity 


C. 6 mil front, infinite back 
shield, 6 mil cell thickness, 
10 ohm-cm base resistivity 

Fig. 5-16 Fluence versus Time 




D. 6 mil front and back shield, 
6 mil cell thickness, 10 
ohm-cm base resistivity 


E. 6 mil front and back shield 
10 mil cell thickness, 

10 ohm-cm base resistivity 


Fig, 5-17 Power Degradation Histories 
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No complete SERT-C cases with shadowing were run, although some shorter 
cases were. A case similar to that listed in the third column of Table 5-4, but 
with the final a = 13200 km (rather than 12800) and i “ 25.3® (rather than 25®), took 
83 days with shadowing and delay (compared to 42 without shadowing). The initial 
orbit has a period of 2. 6 hours. The time in shadow is 36 minutes and the delay 
time is 14 minutes. Together they are almost one third of the orbital period. Thus 
for lower orbits, shadowing and delay in thruster turn-on can considerably lengthen 
the flight time, at higher altitudes proportionally much less time is spent in shadow. 
Since the spacecraft spends more time at lower altitudes, the power degradation 
is more severe, also lengthening flight time. 

Four GEOSEPS-type cases were run. Table 5-5 lists conditions for these 

cases. 


Table 5-5 GEOSEPS Example Data 


Initial a = 

29378 km 

Initial e = 

0 

Initial i = 

7® 

Initial mass = 

2 796 kg 

Maximum initial jet power = 

13 kw 

Specific impulse = 

2900 sec 

Final a = 

42164 km 

F inal e = 

0 

Final i = 

0® 


These cases assumed attitude constraints, power degradation, oblateness, and 
shadowing. A five day time step and a 4-point quadrature between jumps points 
and shadow entrance and thruster turn-on points. Jump points were calculated. 

The launch date was March 21, 5, and two initial nodal angles were considered. For 
each, cases with and without delay were run. Table 5-6 summarizes the results, 

A change in the line of nodes of 45° results in a 6% decrease in transfer time, the 
addition of a delay in thruster turn-on time causes only about a 2% increase in flight 
time. Not all orbits of these transfers enter the shadow and because of the high 
altitude the time in shadow and the delay time is small compared to the orbital 
period (unlike the early orbits of a SERT-C mission). There is little variation in 
for these cases. Because of the initially high orbit, power degradation is less 
than 6%, more than half of which occurs during the first orbit. 
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Table 5-6 GEOSEPS Example Results 


Longitude of the node 

0 

0 

45° 

45° 

Delay 

yes 

no 

yes 

no 

(km/sec) 

.887 

. 885 

. 876 

.874 

Transfer time (days) 

37. 9 

37.3 

35.4 

34. 9 

m./m 
f' 0 

.969 

. 969 

. 970 

.970 

Damage factor 

. 946 

. 947 

. 947 

. 947 

Initial costate 

1179.1 

1161.2 : 

1033.3 

1023.0 


-184. 7 

15. 9 

-101. 1 

-23. 0 


1 -882.2 

-661.4 

-464. 6 

-345.4 


‘ 12643. 

13003. 

-23487. 

-23858. 


-44386. 

-44125. 

-27212. 

-26826. 


-1512. 8 

-1510. 7 

-1292. 1 

-1290. 1 


-1517. 2 

-1533. 8 

-1304. 3 

-1318.4 


The semimajor axis and inclination history are plotted in Fig. 5-18. Their 
changes are monotonic. Fig. 5-19 shows the yaw for the initial and final orbits. 

The yaw for the initial orbit varies fairly smoothly 60^ either side of the orbit plane. 
The spacecraft thrusters are off for eccentric longitudes between 168 and 201° 
because of the shadow. The yaw for the final orbit varies 90° either side of the 
orbit plane and there are two jumps in the yaw angle, one at 151. 6° and one at 
317.4°. 


a 

(10^ km) 



i 

(°) 


Fig. 5-18 Semimajor Axis and Inclination Histories 
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Eccentric Longitude (°) 


Fig. 5-19 Yaw for the Initial and Final Orbits 


The power variation over the initial and final orbits is shown in Fig. 5 20. 
For the initial orbit power reaches . 855 of the maximum with an average of . 938. 
For the final orbit the minimum of . 707 is reached at the jump points. The 
average is .853. 



Eccentric Longitude (°) 


Fig. 5-20 Power Variation for the Initial and Final Orbits 
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The primer vector angle, «, sun angle, g, and thrust angle. 4> are plotted 
for the final orbit in Fig. 5-21, The sun angle varies between 13 and 167° , The 
thrust angle and primer vector angle coincide at 180°. When ct passes through 90° 
there is a jump in ip of 78° when g = 64° (at F = 151, 67° and a jump of 87° when 
|8 = 101° (at F = 317,4°). 
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Fig. 5-21 Primer Vector Angle (a). Sun Angle (/5), 
and Thrust Angie (^) for the Final Orbit 


Fig. 5-22 plots the time in shadow and the corresponding delay time (added 
to shadow time) versus the flight time. Time in shadow reaches a maximum of 
about 59 minutes. Delay is always at least 12 minutes when thrusters shut down. 



Fig. 5-22 Shadow and Delay Times 
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Two other cases were considered that had Lnitlai elliptical orbits. One case 
assumed zero pitch but free roll and terminated at an equatorial geosynchronous 
orbit. The other had both roll and pitch constrained, but was a shorter transfer. 
The spacecraft parameters were the same as for the SERT C cases without 
shadowing discussed previously. Table 5*7 summarizes the initial and final orbits 
and transfer characteristics. 

Table 5-7 Elliptical Orbit Examples 


Roll 

free 

constrained 

Initial a 

10400 km 

10400 

Initial e 

. 378 

. 378 

Initial i 

28. 3® 

28. 3 

Initial Q, 

0® 

0 

Initial (o 

CD 

O 

O 

90 

Final a 

42164 km 

18100 

Final e 

0 

.3 

Final i 

0® 

22. 2 

Final aV 

4.46 km/sec 

1. 72 

Final t^ 

181.3 days 

72.0 

Final m^/m^ 

. 855 

.941 

Final damage factor 

. 627 

. 638 

Initial costate 

9947.3 

9612.4 


4365.8 

7945. 8 


-767. 1 

-470.9 


-718. 2 

-1162. 5 


-16925. 

-19673. 


-26014 

-10090. 


-112. 22 

-77. 15 


Fig. 5-23 shows the semimajor axis, eccentricity, and inclination histories 
for the free roll example. There was also a rotation of q from 0^ to -70® and a 
rotation of to from 90® to 217®, primarily caused by earth oblateness, although 
some perlcenter rotation may have been caused by an effect to avoid intense radia- 
tion, thus lowering power and Increasing transfer time. Power dropped from 
4, 828 kw to 3. 03 kw; most of the drop occurred during the first 70 days as in previous 

examples. 
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One of the interesting characteristics of the pitch and roll constrained example 
was that there were two points per orbit at which the primer vector projection was 
perpendicular to the sun vector (et ~ ±90°) except very close to the end of the 
trajectory where there were four such points. There were no jumps in the thrust 
direction for the orbits with only two such points since the sun angle, g, was not 
between 45° and 135°. The yaw angle is plotted in Fig. 5-25 for the initial and 
final orbits and for a segment of an orbit at 66 days. There are two jumps during 
the final orbit but none for the orbit at 66 days or for earlier orbits. Fig. 5-26 
shows the primer vector angle, a, the sun angle, |g , and the thrust angle, tp, for 
the final orbit. The primer vector angle crosses the 90° point four times, but for 
the first two, g is not in the 45 to 135° range. Then the primer vector angle 
briefly goes above the 90° line between F = 286° and F = 336° . At F = 226° , 
g = 48°, and there is a 30° jump in i/). At F = 336° , g = 85° and there is nearly a 
90° jump in ip. For the 66 day orbit the primer vector curve looks very similar 
but it does not quite reach the 90° line in this orbit segment and these two jumps 
do not occur. 



Fig. 5-25 Yaw for the Initial, 66 Day, and Final Orbits 



M l» tao MC \ 

Eccentric Longitude (^ ) 


Fig. 5-26 Primer Vector Angle (a). Sun Angle ($), 
and Thrust Angle (ip) for the Final Orbit 
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5. 3 Run Time and Accuracy 


Some of the additions to the new version of the program have contributed to 
greater run time, although averaging and the analytic radiation model cause the 
run time to be orders of magnitude shorter than a precision trajectory calculation. 
Many trajectories are run for one optimization, however. Shadow calculations 
cause an increase in run time and the resulting trajectory often requires a shorter 
integration time step to maintain the accuracy needed to obtain convergence. The 
new thruster turn-on delay calculation requires additional time. The constrained 
thrust direction c alculation requires additional calculations which must be performed 
at every quadrature point of every orbit for every trajectory. Calculation of possible 
jumps points in the control greatly increase run time. For accuracy a separate 
quadrature may be used between jump points and shadow entrance and exit points. 

With no shadowing and two jump points there will be two quadrature intervals. 

With shadowing, and if there are two jump points not in the shadow, then there are 
three intervals. There may be up to 6 jump points (for non-circular orbits). (Tump 
points refers to points where « = + 90. There may or may not be an actual jump 
in ^ depending on the value of fl. ) Some cases were run without calculating jump 
points and just arbitrarily dividing up the orbit into two quadrature intervals. The 
loss of accuracy adversely affects convergence and the increase in the number of 
trajectories calculated may increase total run time. For the results in this section 
a four-point quadrature was used. For higher order quadratures run time is 
greater. A smaller time step increases the run time of a trajectory. 

The results in this section were generated on an IBM 360/75 computer. The 
program is coded in Fortran IV on the H compiler. These runs took several minutes 
each. To get a rough idea of the time needed, the approximate CPU time per 
"time step" was calculated. With oblateness, degradation and attitude constraints, 
two 4-point quadratures per orbit and jump points calculated the value was . 025 
min. Thus for a trajectory 250 days long with a 10 day time step there were 25 
steps or .625 minutes per trajectory. Four iterations of the Newton-Raphson 
iterator generate at least 26 trajectories (4 nominals, 3 sensitivity matrices and 
one time history), requiring 16. 25 minutes. Without calculating the jump points 
the CPU/time step was approximately . 018 so the same case would take 11. 7 minutes. 
The shadow cases in this section had similar CPU/time step, often one of the jump 
points was inside the shadow, and so not calculated precisely, reducing computation 
time. 

Experience indicates that using higher order quadratures and smaller time 
steps than used for the results in this section change and transfer time very 
little, perhaps 1 or 2%. However, the initial costate may be changed on the order 
of 10 or 20%. For mission comparisons and to obtain approximate trajectory 
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characteristics less accuracy may be accpetable. If costate histories are to be 
used for guidance, increased accuracy and therefore smaller time steps and 
additional points per orbit must be used. 

More complex and longer cases required better initial guesses to obtain 
convergence. Such cases usually required an iterative procedure involving 
converging to successive points along a reference trajectory. 
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SECTION 6 


CONCLUDING REMARKS 


The program discussed in this report is a generalization of a previously 
developed code which calculated time optimal geocentric trajectories for solar 
electric propelled spacecraft with an optional high thrust stage. The new code continues 
the use of the method of averaging and an analytic radiation and degradation model in 
order to reduce computation time. A nonsingular set of orbital elements is utilized. 
Oblateness, shadowing and solar motion may be included. For the new code, attitude 
constraints may be included, a delay in thruster turn-on upon leaving the shadow is 
modeled, and the radiation and power loss model is revised and generalized to 
various solar array parameters. 

A suboptimal control for the attitude constrained case of zero roll and pitch 
which is nearly time optimal without wasting fuel has been developed. Because of 
the constraints the panels may not always directly face the sun, and power becomes 
a function of spacecraft orientation. The resulting control strategy may call for 
discontinuous jumps in the yaw and panel orientation angles. The minimum possible 
power during an orbit which occurs at a jump point is . 707 of the maximum, although 
the average is usually , 9 or higher. 

The code is computationally fast due to averaging and the analytic radiation 
model compared to precision code (although unlike most precision codes, many 
individual trajectories are run in order to converge to the desired final conditions). 

The code is applicable to general orbits, spacecraft and solar array characteristics. 

The radiation and power loss model represents a close approximation to 
actual data. The analytic format is easy to code and computes rapidly. Partial 
derivatives, necessary for the costate equations, are easy to obtain. The 
coefficients in the mathematical models can easily be altered for new data. Shield 
thickness, solar cell thickness, and cell base resistivity may be varied. 

The time step for the integrator and the number of points calculated per orbit 
for the averaging quadrature affect accuracy and computation time. Experience 
with a limited number of examples indicates that using 8 points per orbit or 16, 
or halving the time step from 10 to 5 days in an over 100 day trajectory may change 
and flight time by 1 or 2%, but the initial values for the costate may change 
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on the order of 10% or more. In order to obtain convergence for complex cases 
better initial guesses and more accuracy is required. It may be necessary to obtain 
a nominal trajectory from a more simple case and then to converge to various 
points along that nominal, at each stage using the initial values from the previous 
stage as initial guesses for the later stage. 

Comparisons were made for a SERT-C type mission between constrained 
and unconstrained cases for an inverse square gravity field. Flight time and 
vary for the constrained cases as a function of launch date and time of day. 

A launch time can be found for which the constrained is less than 1%, and the 
transfer time less than 5% more than the unconstrained case. Jumps in the control 
did occur during the last 20 days of a 130 day mission. For the examples 
considered most orbits had at most two jumps in the thrust direction. Noncircular 
orbits can have up to six jump points. 

The time in shadow and the delay in thruster turn-on can be a significant 
proportion of the orbital period at lower altitudes. As a result the spacecraft 
spends greater amounts of time in the more intense lower radiation zones, thus 
power degrades faster which contributes to longer transfer times. 

The inclusion of power degradation can nearly double transfer time for a 
SERT-C type mission. The intense radiation at lower altitudes causes the power 
to drop to about 50% of its initial value, the exact amount depending on solar cell 
characteristics and the particular trajectory. Most of the degradation is in the 
early portion of the trajectory. 

Further refinements of the code could include extending the valid region of 
the radiation model to those low altitudes between latitudes 50° and 75°, and 
effects to make the program more efficient in terms of individual trajectory 
calculations and convergence characteristics in order to reduce run time further. 

The program offers a tool for mission performance evaluation and compar- 
isons. The costate histories generated are first order approximations to precision 
costate histories and with sufficient accuracy control the code can form the basis 
of a guidance scheme. 
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APPENDIX A 


HIGH THRUST EQUATIONS 


A. 1 Introduction 

In this appendix the details of the high thrust phase are discussed, including 

(12 13 ) 

a brief description of Huntington SmalPs technique and code ' which were used 
for the initial high thrust phase. Since this code uses a special set of variables, the 
transformation to the equinoctial state and costate, used in the following low thrust 
phase, is presented. This appendix is a condensed version of Section 4 of Ref. 2. 

A. 2 The High Thrust Code 

The initial high thrust phase of the program utilizes code developed by 
H. Small, This code calculates time-open minimum trajectories in an inverse 
square field. This subsection discusses the method of Small and presents certain 
equations some of which will be used in Section 4.3 which considers the interface 
between the high and low thrust code. 

The following is taken mainly from References 12 and 13. Taking the orbital 
elements as state variables and the characteristic velocity as independent variable, 
the rates of change of the state variables can be written 

j = 1,..,5 (A. 1) 

in which the unit vector in the direction of thrust, and u, an angular variable 

in the instantaneous orbital plane are controls and the problem is to generate 

extremals which minimize the independent variable between end states. The 
5 

Hamiltonian is H = X d. g and is constant and can be taken equal to unity. The 
j=i 

2 2 

optimal ^ " Ay optimal u maximizes H(u) = [X^(^)] • 
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Using any set of orbital elements 


Av = AiSin(u+A 2 )iR 

+ [a3 + (1 +^())AiCos(u+A2 >] U 1 l 

+ A 4 SLn(u+A 5 )/^ 


(A. 2) 


in which and e. are unit vectors in the radial, circumferential and out- 

of-plane directions; ana = 1 + ecosf where e is the eccentricity and i - u - . 

the true anomaly. The A's are functions of x and X and will be defined later for 
our elements. 

■' 2 

The Hamiltonian will be maximized at some u if [X^lu+Au)] 2 : 0 with 

equality only at ^u = Defining t = tan(Au/2) this inequality 

becomes the ratio of two sixth-degree polynomials in 1 . Small then reduces the 
maximizing test to the form 

t^ j t+a 3^ )> OVt (A.3) 


He then develops a general test for the positive definite quartic. The require- 
ments are {Eq, 4 of Ref, 13) 


ao > 0 

16ci!^(q!2+3z j) > (3 ck^-8o(^ 0!2)^ if 3a^ > 80/^02 

+ 902^2] ' “2^2 “ ^ 


(A. 4a) 
(A. 4b) 

(A, 4c) 


where 

Zl = 4»^a^ - aj «3 

2 2 

Z 2 = Oi«4 + 0io0!3 - ori«2«3 ^ 

There are other simple conditions that the must meet such as 04 > 0 but they 
are all contained in Eq. (A 4). Eq. (A. 4) will be violated first as the coefficients 
change during a firing requiring a switch through an angle described by 

t _ 2(o;^-4a^o;2)(a2'^3zi) - 4 q,^ (9z2-8o(2Zi) 

4«^«3(a^+3z^) + «j(9zj-8«^zj') 

The > 0 case generates an impulse. 
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The H =0 condition (Eq. A.3) and the are explicit functions of eight 
variables in Eq. (A, 2) but because (X^) can be scaled to 1 and u just depends on an 
arbitrary reference axis, It turns out that they can be written as functions of six 
independent parameters. Small chooses e sinf, 0, 0 , T where 

X = sln 0 e + COS 0 cos Te, + cos 0 sLnT e (A. 7) 


(illustrated in Fig. A-1, cos 0 ^ 0) and two others defined by 
k H AiCOs(u+y\2)l ) 

j = |a 4 s^*^(^"^A 5 )sin 0 + A 4 ^‘^s(u+A 5 )cos 0 cosT 

+ [AjCOs(u+A2^® 

« A cosf^ cos 0 sLnT| /cos ^0 

= |^A 4 COs(u+A 5 )/<^^S 5 Z^ esin fsinT^cosT 

+ [tan 0 - ke sin f]sin T 
Further manipulations yield 


tan 0 


(j sin T -i- ke sinf) 
1 + k 0 cos T 


and expressions for the 


“o 

= ♦[: 

1 - - 

(2+^k^)tan^<z5 J - 

“1 

= 4ji) 

tan 0 ( 2 k- 

2 

COST )- 2e sinf(l-2 tan 0) 

“2 

" “o 

+ «4 + 8 

esin f tan 0 (cos T-k) 



- 4(0-l)| 

cos T-k)^ - tan ^0 J 

“3 

= «1 

+ -Dtan 0 (cos T-k) 



+ 4 e sinf 

j^(cos T-k)^ - tan^ 0 j 


«4 = [(3-!/))cosT - 2k][2k-cosT] - ( 0 -l)sin^T 


(A. 8 ) 


(A. 9) 


(A. 10) 


(A. 11) 
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Figure A-1. Firing angles T and 0- 

The k and j variables can also be written In terms of the radius and 
velocity vectors and their adjoints. 


k = * R)/h cos 0 

(A. 12) 

ah \ 

4 - V 

J ' 7^ 

h cos ^ 

(A. 13) 

where h is the angular momentum magnitude and h is 
magnitude at some reference orbit, u is the gravitation 

the angular momentum 
constant, and 


(A. 14) 


which is a constant of motion. 
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Small notes some special cases* If e = 0, ‘ <^3 and 01 2 ~ ^ 4 

that = 0 and Eq. (A.. 4c) reduces to > 0. Eq* (A. 4a) is automatically satisfied. 

The double maximum occurs when z = 0 at t ^ = - 2ry /ao. A second case Is the 

1 s o' 6 

Generalized Hohmann Transfer for which esLnf = j 0 Implies tan0 “ “ ^3 " 

Z2 = 0 so Eq. (A. 4c) reduces to > 0 and Eq. (4.4b) to ^ ^2 This 

equation turns out not to be restrictive and the double maximum occurs when 
= 0 at ^u = 180° . 

Small's code works In the following manner. An initial set of parameters 

-1/2 1/2 

(e, f, ) which satisfy Eq. (A 4) is input (for our program e = 0, 

f = 0, ^ = 1 so the input is (0, 0, T, j, k)). The initial h is arbitrary, just 

scaling t = h j \x which is the independent variable within the code. We have set 

h- = h^ = 1. Then t is incremented by a succession of The corresponding 

values of , 0 and Eq. (A. 4c) are computed until Eq. (A 4c) is zero at some The 

values of AiSin(u + A ^cosiu-^/^ A 4 sln(u+A3) and A4 cos(u-f A5) are 

computed from the known values of e, f, T , j, k at t and all variables are 

s 

switched to a new peak by adding ^u computed from Eq. (A. 6); then the next impulse 

s 

continues by incrementing t again. The result is a series of known Impulses and 
coasting arcs. Since for our case the initial orbit was circular, the location of 
the first impulse was arbitrary (as far as the high thrust code is concerned) and 
so the initial f = 0. We limited the number of impulses to two and specified 
the total ^V. In this case only the remaining available was used on the final 
impulse. The location was optimal although the might be less (or more) than 
called for by Small's program. 

A. 3 Coupling of High and Low Thrust Code 

In the combined high and low thrust program the output of Small's code must be 
interfaced with the beginning of the low thrust code. In addition we have constrained 
Small's code by 1) requiring the initial (high thrust) orbit to be circular, 2) 
specifying the maximum number of initial impulses and the corresponding total ^V. 

The first constraint allows us to calculate the limits on the permissible input to 
Small's code, i.e., the requirement that the input satisfy Eq. (A 4) 

The requirement that the initial high thrust orbit be circular and the 
assumption that the first impulse is in the direction In order to raise the orbit, 
combined with the inequalities of Eq. (A. 4) yield the following requirements on T, 
k and j. 

I T I S I (A. 15) 
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(A. 16) 


— - < k < cos T 


Ij! < (l+kcosT) J (A. 17) 

cos T 4 k 

These constraints guarantee that values picked for T ^ k, j will yield an optimal 
impulse. 

For a given total aV and a given set of initial conditions, a one, two, three 
or more impulses optimal transfer may result. We are limiting the number of 
impulses to two. If the specified total aV is less than the maximum for the first 
impulse calculated by the code, that impulse is still optimal under the constraint. 

If the total aV is greater than the maximum for the first impulse, then all of the 
remaining aV will be used at the second impulse. If the remainder is less than or 
equal to the maximum for the second impulse calculated by Small^s code, then the 
transfer is still optimal under the constraint. If the remainder is greater, then 
the transfer is not optimal. A three or more impulse transfer could have reached 
the same orbit with less Usually for transfers of interest, one or two 

impulses will be optimal. 

For the calculations in Small* s code a *'first orbit" coordinate frame was 
assumed. In this frame the initial orbit has zero inclination and the initial impulse 
is assumed to occur at f = 0. Also o) ' O “ ^ initially. The orbital elements and 
their adjoints which are contained in SmalUs S(I, J,K) array are with respect to 
this coordinate system. The elements of the S array are made up of orbital 
elements and other of Small's variables. This array is calculated by Small's code 
and it is from this array that we must obtain the equinoctial orbital elements and 
their adjoints and other information for input to the low thrust code. In Small's 
code the 1 in S(I, J, K) varies over different trajectories (he calculated neighboring 
trajectories along with a nominal trajectory for a purely high thrust Newton 
iteration). For our program I was set to 1. J varies over different orbits. For 
our case J varies from 1 to 3 for a two impulse trajectory. K varies from 1 to 20 
and these quantities are listed in Table a. 1 In what follows, for brevity, 
will be used for S(I, J, K). 
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Table A-1+ 


K 

S(I, J,K) 

11 

j/)^/^(cosT 

1 

T 

12 

COS(u + Q ) 

2- 

DT 

13 

sin(u + n ) 

3- 

Ip ^^^tan(^u/2) 

14- 

cos 

4 

(h/h-))/j ^'^^cos0 

15=c 

sin ^u 

5 

-1/2 . , 
Ip ' esinf 

16 

cos L 

6 

r , 1/2 

k ip ' 

17 

sin i cos u 

7 

tan 0 

18 

sin i sinu 

8 

^4COs(u+A5)/cos 0 

19 

sin i cos 0 

9 

i^l/^sinT 

20 

s in i s in f). 

10 




Information for impulse 

is stored in S(I, J+1, K). 


+This 

is essentially Table D. 

1 of Reference 12. 


In our 

program, h'’' was set to 

one and the initial orbit 

is assumed 


Thus for each orbit (J = l, 2, 3) 

h = /cos 0 


u2//i 2 . 

a - a^h /(1-e ) 


c = ''’^10 ‘^ 5 ]^ 

In the "first orbit" coordinate frame, 

cos i = 

1 


sin L 


^ 


, „ - 20 
tan Q ~ "C — 

^19 


In addition 


tan u = 


tan f = 


®ia/®i7 

Sio-1 


- k) 


circular. 

(A. 18) 
(A. 19) 
(A. 20) 

(A. 21) 
(A. 22) 

{A. 23) 

(A. 24) 
(A. 25) 


102 



05 

= u - f 

(A. 26) 

tan 0 

S 7 

(A. 27) 


Sn 


tan T 

9 

(A. 28) 





Now the orbital elements of the first orbit are known with respect to the 
inertial reference coordinate frame, a, e, ij, f)j, (jj. (obtaining is discussed 
later). In the "first orbit" frame, i^ = 0, Q ^ = 0, = 0. A vector in the first 

orbit frame is related to a vector in the reference frame by the transformation 





COSic^COSQ^ I 

COS tjo^sinOj 

sin ^ 3 ^ cos ij 

- sin cos ij sinp^ 1 

+ sin cos vOj 


- sin Oj cos 0^ 1 

- sin ojj sinOj 

cos sin i^ 

- cos j cos ij sin Oj j 

+ cos 00 j cos ij cos p j 


sinijSinPj 

- sin ij cos p ^ 

cos ij 

_ 




(A. 29) 

The coordinate frame dependent orbital elements of the second or third 
orbit are obtained in the "first orbit" frame by using the S array; call these 

i . D » coo- The orbital elements in the reference frame, ig, p.g, 0 ) 3 . can be 
2 2 2 
obtained from 


Tg = TgTj 


(A. 30) 


where Tg is given by Eq. (A. 29) with the subscript 1 replaced by 2. The elements of 

T„ are defined similarly. T„ can then be solved for ig, Qg, 0 ) 3 * Let be the 

f h ^ 

L - j element of Then 


sin 



32 


(A. 31) 


cos ig /^gg 


(A. 32) 
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(A. 33) 


sin 103 " - 1 ^ 3/310 ig 
cos (Og ' ■^• 23 /®'^^ 4 (A. 34) 

sinpg = t 3 j/sin ig ^A.35) 

cosQg = - 4 (A, 36) 


If = 0 then a ;3 can be found from 



cos((j^j 2 ” "^11 

(A. 37) 


sin(^2 ^3^ ~ ^12 

(A. 38) 

The equinoctial elements, h. 

k are given by 



h = e sin(o )2 O 3 ) 

(A. 39) 


k = e cos((jt )3 + 03 ) 

(A. 40) 

It is possible to get p, q directly from T^. 

^31 

(A. 41) 


p = 

^ ■*' -^33 


, . ■'•32 

^ ■^33 

(A. 42) 

The A'*® given in terms of the classical state and costate by Eq, A -8 

of Ref. 12. This can be inverted to give the costate in terms of the 



[a3 ' e AiCos(A2+o,)] 

(A, 43) 

K = V 


(A, 44) 

>-i = 

’ n 2 , A4SinA5 
ad-e^) 

(A. 45) 
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^ + to) 

w ^a(l-e2) ^ ^ 


(A. 46) 


X 

a. 


/— H 

^ aU-e^) 


^ .cos^ j-sin L + X cos i 

4 0 CO 


(A, 47) 


However, this costate is given in terms of the first orbit coordinate frame. Since 
in this frame i = 0 and e = 0 initially, then initially X = X = 0. In the equatorial 

CO 44 

reference frame if Q and co ^.re free, then = 0 and X^ - 0 initially. The X^ 

and X will be the same in the equatorial reference frame as in the "first orbit 
e 

frame". The form of X. at the Lnltial time follows from the requirement that X =0. 

L 44 

This condition also yields the location of the first impulse. If in the reference frame 
fj is fixed initially then x„ is free though \ , is still zero. 

In the reference frame the first impulse occurs at The adjoints to p 

and i in the reference frame can be realted to {i)j by the following equations derived 
from results listed in Appendix B of Reference 2. 


X/ 



(A. 48) 


A4SinijCOs(A^ - Wj) 

^o 


(A. 49) 


If in the reference frame, we require that X' = 0, then for sin i 7 ^ 0 


A5 - 0)i 


- ± 


1 

2 


or 




^5 
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Taking the convention that X/ and have the same sign determines that 


Then 


It 


(A. 50) 



(A. 51) 


If the initial q is fixed then X is free and not necessarily zero although 
is still zero. For this case it is no longer necessarily true that ” A5 "“j * 
rather ^oj arbitrary but then determines X^ and \/ through Eqs. (A. 48) and 
(A. 49). 


The formulas for X and X are valid at any time during the initial high 
a e 

thrust phase. During the high thrust phase X^ remains constant, but X. and X 

Q I to 

may change. Initially, there are some ij, flj, then using results from 
Appendix B of Ref. 2 and the fact that X' =0, let 


Cl = X/ cosQj - X^cotiiSinQi 

Cg = X/ sinfii + X^cot ijcos Qi (A. 52) 


For any set of i, (j, to tn the reference frame 

= Cl cos n + Cg sin Q 

n 3 

X = sin i(c- sin Q - c« cos O ) + c^ cos i 

tli I ^ o 


(A. 53) 
(A. 54) 
(A. 55) 


For SEP and high thrust in an axially symmetric field, the fixed q option 
should be used (in which case the initial ui should be driven to ^ so as to 

make = 0). If this option is not used and o is allowed to be free initially, then 
varying fi initially will have no effect on the final conditions, and the partial deriva- 
tive matrix will be singular. 

We need ^ and in terms of S array. From Eq. (4.2), (4. 7), 

(4. 8) and from Table B-1 the following quantities are obtained. 
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^3 

COS )k - cos T ] 

(A. 56) 

A1COSA2 = 

cos 0[cos T-k]cosu + sin 0 slnu 

(A. 57) 

AfSinAg = 

-cos 0[cos T“k]sinu + sin 0cosu 

(A. 58) 

A 4 COSA 5 = 

cos 0 (Sg cos u + sinT sin u ) 

(A. 59) 

A4sinA5 = 

cos 0(-SgSinu + ip sinT cosu) 

(A. 60) 


In order to calculate X and X we need the combination 

a e 

/ ^11 \ 

cos cos to - Aj sin ^3 sin to = sin f sin 0 + cos f cos 0 ( ■ ■ ) (A. 61) 


where we have used u = o) where ^ and f were derived in terms of the S 

array previously. Obtaining A3 terms of the S array, 

A 3 = T''4®10-S6-Su’ 

is needed only on the first orbit where u = f - ^3 ' on the first orbit 

A^cosA^ ~ Sg cos 0 (A. 63) 

A4sinA5 = Sg/S^cos0 (A. 64) 

A 4 can be obtained by taking the square root of the sum of the squares of the 
above two quantities. There is an ambiguity of sign, however. The computer code 
picks the sign by the criterion that if the inclination is desired to be reduced, Aq 
(and thus X ) should be negative; similarly if it is desired to increase inclination, 

X. and A4 should be positive. This choice can be overriden by an input to the 
program. Once the sign of A4 is picked, A5 can be found unambiguously (as long 
as A4 t ^ i calculated. 
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Finally the equinoctial costate is given, at each orbit using Eq, (3. 25), 




(A. 65) 


sin + cos {^^+ q ^) ^sL 

(A. 66) 

\ = 

cos(t03+O3)X^ - 310(103+03)^^ 

(A. 67) 

^ = 

2 ^3 cosOo 

2sinn cos 7 r^i+ ^(X_-X ) 

tan(l2/2) ^ ^ 

(A. 68) 

X = 
q 

i-o s In rt 

2cos0 3Cos‘^-4.X - ^ ) 

2 1 tani /2 ^ w 

(A. 69) 


A multiplicative scale factor may also be needed before input to the low thrust code. 

Since and therefore cannot be calculated unambiguously if =0, 
this case can lead to difficulties in running the program. “0, 

This implies a zero inclination change during the initial high thrust phase. Because 
of the radiation and shadowing effects, which for SEP will generally cause an 
optimal trajectory to contain some inclination variation even if the initial and final 
inclinations are the same, this case will generally not arise for SEP. It could 
arise for constant power and zero inclinations change and so the program may not 
work for this case. The initial guesses for T and j should not both be zero. 
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APPENDIX B 


THE MATRIX M AND ITS PARTIALS 


Table B -1 Elements of M 


2 X, 


2 Y, 


M 


11 


M 


21 


M 


22 


M 


23 


M 


31 


M 


32 


M 


33 


- -j— , jvi^2 

- , IVl. 

n a 

n a 


ax X, 

2 

na 

ak n 



2 


na 

bk n 

k(qYj-pX^) 


2/1 u2 
na yi-h -k 


ax^ Xj 

2 

na 

bh n 


"aY^ 

2 

na 

n 

- h(qyj-pXj 

) 

2 L ,2 ,'2 

- 

na yi-h -k 



13 


= O 


(sinF-h B ) 


M41 O, M42 " ^43 


(l+p^+q^)Yj 


2naVl-h^-k^ 


M 


51 


o. M52 = o. M53 


(l+p^+q^)Xj 


2na^yi-h^-V 


109 



Table B-2 PartLals of and with 
Respect to h and k 


hh 

Bh 

BXj 

Bk 

BY 


1 


Bk 


2 3 

-hpcosF~(g+ h ) (hcosF-ksLnF) 
kgcosF-1 + hkg^ (hcosF-ksinF) 
hgsinF-1 - hkg^ (hcosF-ksinF) 
-kgsinF+(g+k^_£^) (hcosF-ksinF) 


Table B-3 Partial of M with Respect to a 


3M _ 1 

aa 2a 


3 0 0 0 0 
0 10 0 0 
0 0 10 0 
0 0 0 1 0 
0 0 0 0 1 


M 


no 



Table B-4 Partial of M with Respect to h 


ah n^a ah ah n^a 9h Bh 


O 


aM 


21 


ah 

BM 


-hMai 

7TT7^ 

l-h -k 


JT:^ 


na 




22 ^ 


-hM 


22 




dh 


; j~7i 

l-h “k 


na 




3^Y, 


^hBk 


1 ^ ^^1 (sLnF-h)9) _ ^ (g+h^^^) 


Bh 

B^ 


n 


1-/9 


2_3, 


Bh 


1 (slnF-hg) _ _1^ (g+h g ) 


n 


1-^ 


BM23 

hM23 

4 k - 

- u 

BYj _ 



Bh 


naVl-h^-' 

-J V 

< 

bh 

Sh / 


BM 31 

-hMsi 

yi-h^-k^ 

B^Xj 

Bl^j 

(cosF-kg) _^^1 

hkp^ 

Bh 

■ , .'S'—S- 
l-h -k 

2 

na 

Bh^ 

Bh 

n n 

1 -P 

aMs 2 

-hM32 



BY, 

(cosF-k/9) ^ 1 

hkg^ 

bh 

1 .2 .2 

l-h -k 

2 

na 

Bh^ 

Bh 

n n 

_ 

BM 33 

hMss 

. ^33 

h 

- fci 

BY, BX, \ 

- p 1 


Bh 

■; ,2 “2 
l-h -k 

h 2 A . 2 . 

na yi-h -k 

? r 

Bh Bh ' 



BM 


41 


BM 


O. 


42 


Bh 


Bh 


BM^3 _ hM^3 
= O. - 2T^ 


Bh 


l-h‘'-k‘ 


Y, Bh 


^^51 

Bh 


BM 


O . 


52 


Bh 


O 


^^53 

Bh 


™53 

^ + 

l-h^-k"^ dh 


111 



Table B-5 Partial of M with Respect to k 
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Table B-6 Non-zero PartLals of M with Respect to p 
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Table B-7 Non-zero Partlals of M with Respect to q 
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Table B-8 Partials of and with Respect to h and k 
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APPENDIX C 


AN ALTERNATE FORM OF THE M-MATRIX 


bz 

In Reference 28, page 30, a new form of the M or — p- equation was written 

with respect to the equinoctial orbit frame. The partials of a, h, k, p and q are 
repeated here. The "M-matrlx" is just made up of the coefficients of the basis 
vectors in the following expressions. 


c»r 


T 
n a 


^ = - [Gf + rX y] +;^ [qY - pX ]w 
^r u " G — 


^ = ^[Gi+rYi]-^[qY -pXj]w 


(C. 1) 


SP = (1+p^+q^) Y ^ 

9r 2G 1 - 


^ w 

9£ 2G 1 - 

where ^ w are the equinoctial basis vectors defined In Section 3. 2 and X^ and Yj 
are the f and g components of position respectively. Also 


G = na 


7 n n- 

1-h -k = YjXj - XjYj 

(C. 2) 

xr/r = (Xj|-Yjf)/r 

(C.3) 

- - 

(C.4) 


(C. 5) 
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A useful form of — for the constrained case when thrust must be orthogonal 
to the radial direction is Eo write it with respect to the e , e , e frame where 

— _ g — w 

e^ is parallel to the radius vector, is normal to the orbit frame and e^ completes 
the orthogonal triad. In terms of the equinoctial basis vectors 

ir = (C.6) 

ig = (- f + Xj g)/r (C. 7) 

= w (C.8) 


Applying the inverse of this transformation to the previous ^ equations yields 
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\ 1 r J—s J 

GYj 

^ / 

, . GX, ^ 
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~r + (rY. + ^)e 1-^(qY. -pXJe 

.. L r -r \ 1 r y — sj G ^ 1 ^ 1 --w 




2G 


^li 


w 




1+p^+q^ 

“TO 


X.e 
1— w 


(C.9) 


When thrust is constrained to be orthogonal to the radius vector, the e 

— r 

terms drop out. If this form is used in evaluating state equations there would be 
fewer calculations. The partials needed for the costate equations would also be 
simpler. The unconstrained case would still involve calculating all the terms. 
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APPENDIX D 


SUN-EARTH RELATION 


For the solar electric program, the sun*s location in equinoctial coordinates 
is needed. The distance from earth to sun varies due to the elllptlcity of earth's 
orbit. The following is taken from Reference 2. 

The sun’s location in earth equatorial coordinates Is given by 

- cos (co+ v) 

rJ = - cos O sin (co +v) 

« sin O sin (co + v) 

where 

O is the obliquity of the ecliptic, 

60 Is the argument of perihelion, 

V is the true anomaly of the earth at 
the time of interest. 

The true anomaly can be approximated by (see Ref. 37 , p, 55) 

3 

V = M + (2e - — )sin M + -e^sln 2M + 3M (D. 2) 

4 4 12 

and M = nt + 

The orbital elements of the earth are taken from Ref. 37, p. 378, for an epoch 
of Julian Date 2436935.0 (1960 Jan. I.5.E.T.). In particular, 

60 = 102?25253 

M = - 1?90562 
o 

e = .016726 
n = .985609^ /day 
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The distance to the sun (in A.U. ) is just 




1 - e 
1+ecos V 


(D.3) 


The difference between the sun-spacecraft distance and the earth-sun distance 
is assumed negligible. 

In equinoctial coordinates 


R 
— s 


[f g Rg 


(D.4) 


where the equinoctial basis vectors were given in Eq. (3.4). In the analysis, it is 

necessary to have the partials of the £ and g components of A with respect to the 

equinoctial orbital elements. These two components are X = R . f and 
^ s — s — 

Y = R . g . Since f and g are functions only of p and q, we need 


ap 

ax 


s 


aq 

aY. 
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SP 

aq 


= R„ 


^s 
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ap 
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aq l+p“+q 

a£ 2qX 

ap 

A 

ag 
aq 


(qY +Z ) 

l+p"'+q'‘ ® ® 


2PYg 

Y7T 


s 

2._2 


(D. 5) 


1+p +g 


which utilize the formulas 


af_ 
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l+p^+q 
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II 

2p 
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.".2,2 
1+p +q 
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1+p^+q^ 

II 

2 

sq 
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APPENDIX E 


WORST CASE EVALUATION OF 
SUBOPTIMAL APPROXIMATION 


In Section 3.6 two reasons were given tor utilizing a suboptimal approxima- 
tion where only part of the Hamiltonian is maximized with respect to the thrust 
direction. The first reason is that it reduces the optimality condition from the 
root of a 6th to a 3rd order polynomial. The latter can be solved in closed form 
with consequent savings in machine time. The second reason is that it is aesthet- 
ically unsatisfying to bias the control so that the desired acceleration component is 
decreased while the thrust is increased. This result is mathematically time 
optimal because it increases the mass flow rate and reduces the mass. While time 
optimal solutions are of more practical interest than fuel optimal solution, it does 
not seem desirable to simply throw fuel awa3^ 

This appendix is intended to show that the suboptimal approximation is 
negligibly different from the time optimal case for practical missions by analyzing 
the worst possible case. This case occurs when the radius vector to the spacecraft 
is always at right angles to the Earth-sun line and the primer vector is orthogonal 
to the Earth-sun line. A simple example of such a case is the coplanar enlarge- 
ment of a circular orbit normal to the earth sun line. Such a case will be 
characterized by a 41% increase in and a 100% increase in flight time as 
compared to the unconstrained case. For the more practical examples treated in 
the body of this report, the fuel and time penalties due to attitude constraints are 
much smaller than this. For these practical cases, the difference between the 
suboptimal and minimum time solutions will be smaller than for the extreme case 
treated herein. 

This case may be characterized by only two state variables, the mean orbital 
speed V and the mass m. The rates of change of these two variables are given 
in Eqs. E. 1 and E. 2. 


2P 


V = 


me 


cosil) sini/j 


(E. 1) 
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2P 

m = - — 


(E.2) 


The Hamiltonian of this problem is given by Eq. E, 3, 


H = - X. 


2P 


me 


cos^fj sinit) - X. 


2P 


m 


COS0 


(E.3) 


The value of the yaw angle «/; that maximizes H is given by Eq, E. 4. 


SLn 


^ A \ ^ ^ A ^ ^ 


4X c 


4\ c 

A. V 


(E.4) 


The rates of change of the Lagrange multipliers are given by the canonical 
Eqs. E. 5 and E. 6. 


^ = 0 


(E. 5) 


2P 


V " ■ V ~2^ 


(E. 6) 


m c 


Three first integrals result from equations E, 5, E. 7 and E, 8. Equation E. 7 
follows from the autonomous nature of the problem while Eq. E. 8 follows from 
Eq. E.3. 


H = 0 


= mX + mX 

m m 


H 


dt 


(E. 7) 
(E. 8) 


The three first integrals resulting from these equations are given by E. 9, E. 10 
and E, 11. The values of the integrals follow from the transversal! ty conditions 
for a minimum time problem. 


H = - 1 


X m = - (t-t.) 
m f 


(E. 9) 
(E. 10) 


IIljC 


(E. 11) 
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Substituting into Eq. E,4 yields E. 11 


s Lnjf) 


t-t^ 


4m^c 


t-t. 


o/ 2 . ^ o 


(E. 12) 


4m^c 


The problem can now be solved by quadrature to yield Eqs. E. 13, E. 14 and E. 15. 


m 

m. 


2cos \j) cot^/^ 


(E. 13) 


Vf ' ^ 


= cr/T-2sm*t3tn52£iLL!S!!iL] 

L-^ 1 + J 


(E, 14) 


me 

t^ - t = [4sLnJ^) - 2CSC0] 


(E. 15) 


The corresponding results for the suboptimal approximation used in the paper are 
given by Eqs, E. 16 and E. 17. 


c m 

V - V = — — 

^/2 


2 

t - t = — (m - m-) 


(E. 16) 


(E. 17) 


The optimal and suboptimal cases are compared in Fig. E-1, which plots 
the normalized time and mass against the normalized velocity change. The plot 
shows that the difference between the two cases is negligible for mass ratios 
characteristic of SERT-C (> .84). For smaller mass ratios the suboptimal 
approximation will take longer but use less fuel than the minimum time case. For 
practical missions the differences should be indistinguishable. 
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Fig* E-1 Normalized Time and Mass Change Versus 
Normalized Velocity Change 
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APPENDIX F 


SINGLE AVERAGED OBLATENESS EQUATIONS 


In these equations is the earth’s radius, 4 is the gravitational constant, n 
is the orbital angular speed, the oblateness coefficient and a, h, k, p, q the 
equinoctial orbital elements. 


Table F -1 J2 Variation of Parameters Equations 


3uRgJ2k[l-6(l^+q^) + 3(p^+q^)^] 
2na^(l-h^-k^) (1+p^+q^)^ 

3pRgJ2h[l-6(p^+q^) + 3(p^+q^)^] 
2na^(l-h^-k^)^ (1+p^+q^)^ 


3 pR^J 2 q(lV-q^) 

2 na^(l-h^-k^)^ U+p^+q^) 


■q- 


SfiR^jgPdV-q^) 

2na^(l-h^-k^)^ (l+p^+q2) 


Table F -2 Partial of J2 Equations with Respect to a 

9h _ . Z — 

ia 2 a 

^ = . I ^ 

Ba 2 a 

^ = - I £ 

Ba 2 a 

ii = - Z i 

Ba 2 a 
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Table F-3 


Table F-4 


Table F-5 
BP 
Bp 
BP 
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of Jg Equations 
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li 
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Partial of Jg Equations with Respect to k 
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Partial of Equations with Respect to p 
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APPENDIX G 


SHADOW CALCULATIONS 


This appendix contains a summary of results from Ref. 2 for shadowing. 
From geometrical considerations an equation can be derived which the entry and 
exit angles must satisfy. Such an equation Is given In Reference 38 and the equation 
given In this section is essentially the same, except that It Is given in terms of 
equinoctial orbital elements. 

The spacecraft position is given by 

where Xj and Yj were given in Eq. (3. 19) and (3.20). Let the unit vector from the 
earth to the sun be given by 


R 

— s 



This is in terms of the equinoctial coordinate frame and thus depends on the 
equinoctial orbital elements p and q. The calculation of the sun's direction in the 
equinoctial coordinate system is discussed in Appendix D. If a designates the 
earth's radius, the cosine of the angle between r and R^ is given by 

^ ^ <l£l^-a2)l/2 (G_2) 


or. 


+ (G.3) 

Squaring and rearranging 

S = (l-x2)x2 . (1 -y^)y2 - . 0 (G.4) 

This is the shadow equation which must be satisfied by the entry and exit angles, 

Xj and Yj are functions of cosF, slnF, a, h, and k (see Eq. (3. 19) and (3. 20)). By 
further manipulations one can derive a quartic equation in cosF. The coefficients 
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of this quartic equation are given in Table G-1. Spurious roots can be eliminated 
by the criteria that S = 0 and that H • £ < 0. In addition, for the entry angle 
^S/bF < 0 and for the exit angle > 0. 

Table G-1 The Shadow Quartic Equation 
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Derivatives of F and S 


The derivative of F with respect to z is needed to evaluate the costate 
equation. It can be obtained implicitly from the shadow equation, 

^ /9S 

dz ^F 

These partials are listed in Table G-2. Note that in calculating ^S/^p and ^S/^q 
we have taken into account the fact that the sun's direction is given in equinoctial 
coordinates and therefore made use of the partial given in Appendix D. 

Table G-2 Partials of the Shadow Function 
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